Index: trunk/Mars/msimreflector/MFresnelLens.cc
===================================================================
--- trunk/Mars/msimreflector/MFresnelLens.cc	(revision 19593)
+++ trunk/Mars/msimreflector/MFresnelLens.cc	(revision 19593)
@@ -0,0 +1,273 @@
+/* ======================================================================== *\
+!
+! *
+! * This file is part of CheObs, the Modular Analysis and Reconstruction
+! * Software. It is distributed to you in the hope that it can be a useful
+! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
+! * It is distributed WITHOUT ANY WARRANTY.
+! *
+! * Permission to use, copy, modify and distribute this software and its
+! * documentation for any purpose is hereby granted without fee,
+! * provided that the above copyright notice appears in all copies and
+! * that both that copyright notice and this permission notice appear
+! * in supporting documentation. It is provided "as is" without express
+! * or implied warranty.
+! *
+!
+!
+!   Author(s): Thomas Bretz,  6/2019 <mailto:tbretz@physik.rwth-aachen.de>
+!
+!   Copyright: CheObs Software Development, 2000-2009
+!
+!
+\* ======================================================================== */
+
+//////////////////////////////////////////////////////////////////////////////
+//
+//  MLens
+//
+//////////////////////////////////////////////////////////////////////////////
+#include "MLens.h"
+
+#include <fstream>
+#include <errno.h>
+
+#include "MQuaternion.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+ClassImp(MLens);
+
+using namespace std;
+
+// --------------------------------------------------------------------------
+//
+// Default constructor
+//
+MLens::MLens(const char *name, const char *title)
+{
+    fName  = name  ? name  : "MLens";
+    fTitle = title ? title : "Parameter container storing a collection of several mirrors (reflector)";
+
+    fMaxR = 25;
+}
+
+// --------------------------------------------------------------------------
+//
+// Return the total Area of all mirrors. Note, that it is recalculated
+// with any call.
+//
+Double_t MLens::GetA() const
+{
+    return fMaxR*fMaxR*TMath::Pi();
+}
+
+// --------------------------------------------------------------------------
+//
+// Check with a rough estimate whether a photon can hit the reflector.
+//
+Bool_t MLens::CanHit(const MQuaternion &p) const
+{
+    // p is given in the reflectory coordinate frame. This is meant as a
+    // fast check without lengthy calculations to omit all photons which
+    // cannot hit the reflector at all
+    return p.R2()<fMaxR*fMaxR;
+}
+
+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;
+    }
+};
+
+float 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) * 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)   ( 0 )
+
+   // 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)
+        return 0;
+
+    const double sqrtBac = sqrt(B2 - a*c);
+
+    if (a==0)
+        return 0;
+
+    const double dz[2] =
+    {
+        (sqrtBac-B) / a,
+        (sqrtBac+B) / a
+    };
+
+    // 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;
+}
+
+double CalcRefractiveIndex(const double &lambda)
+{
+    // https://refractiveindex.info/?shelf=organic&book=poly(methyl_methacrylate)&page=Szczurowski
+
+    // n^2-1=\frac{0.99654λ^2}{λ^2-0.00787}+\frac{0.18964λ^2}{λ^2-0.02191}+\frac{0.00411λ^2}{λ^2-3.85727}
+
+    const double l2 = lambda*lambda;
+
+    const double c0 = 0.99654/(1-0.00787/l2);
+    const double c1 = 0.18964/(1-0.02191/l2);
+    const double c2 = 0.00411/(1-3.85727/l2);
+
+    return sqrt(1+c0+c1+c2);
+}
+
+void ApplyRefraction(MQuaternion &u, const TVector3 &n, const double &n1, const double &n2)
+{
+    // u: incoming direction
+    // n: normal vector of surface
+    // n2: refractive index of new(?) medium
+    // n1: refractive index of old(?) medium
+
+    const double r = n2/n1;
+    const double c = n*u.fVectorPart;
+
+    const double R = 1 - r*r*(1-c*c);
+
+    const auto v = r*c + sqrt(R<0 ? 0 : R);
+
+    u = r*u - n*v;
+}
+
+Int_t MLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const
+{
+    const double F = 50;    // [mm] focal length
+    const double R = 25;    // [mm] radius of lens
+    const int    N = 25;    // [cnt] Nuber of grooves within radius
+
+    const double w = R/N;   // [mm] Width of a single groove
+
+    const double n = CalcRefractiveIndex(wavelength==0 ? 450 : wavelength);
+
+    const int nth = TMath::FloorNint(p.R()/w); // Ray will hit the nth groove
+
+    const double r0  =  nth*w;      // Lower radius of nth groove
+    const double r1  = (nth+1)*w;   // Upper radius of nth groove
+    const double cen = (nth+0.5)*w; // Center of nth groove
+
+    // FIXME: Do we have to check the nth-1 and nth+1 groove as well?
+
+    // Angle to the center of the groove
+    const double phi = acos(cen / (F/1.125));
+
+    //        theta                peak_z
+    const Groove g( TMath::Pi()/2 - phi, r0*tan(phi));
+
+    // Calculate the insersection between the ray and the cone
+    float dz = CalcIntersection(p, u, g, -3);
+
+    // No groove was hit
+    if (dz>=0)
+        return -1;
+
+    // Propagate to the hit point
+    p.PropagateZ(u, dz);
+
+    // Check where the ray has hit
+    const double pr = p.R();
+
+    // If the ray has not hit within the right radius.. reject
+    if (pr<=r0 || pr>r1)
+        return -1;
+
+    // Normal vector of the surface at the hit point
+    // FIXME: Surface roughness?
+    TVector3 norm;
+    norm.SetMagThetaPhi(1, g.fTheta+TMath::Pi()/2, p.XYvector().Phi());
+
+    // Apply refraction at lens entrance (change directional vector)
+    // FIXME: Surace roughness
+    ApplyRefraction(u, norm, 1, n);
+
+    // Propagate ray until bottom of lens
+    const double v = u.fRealPart; // c changes in the medium 
+    u.fRealPart = n*v;
+    p.PropagateZ(u, -3);
+    u.fRealPart = v;  // Set back to c
+
+    // New normal surface (bottom of lens)
+    norm.SetMagThetaPhi(1, 0, 0);
+
+    // Apply refraction at lens exit (change directional vector)
+    // FIXME: Surace roughness
+    ApplyRefraction(u, norm, n, 1);
+
+    // 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());
+
+    // This is only correct if the focal distance is calculated from the inner lens surface!
+    p.fVectorPart.SetZ(0);
+
+    return nth; // Keep track of groove index
+}
Index: trunk/Mars/msimreflector/MFresnelLens.h
===================================================================
--- trunk/Mars/msimreflector/MFresnelLens.h	(revision 19593)
+++ trunk/Mars/msimreflector/MFresnelLens.h	(revision 19593)
@@ -0,0 +1,32 @@
+#ifndef MARS_MLens
+#define MARS_MLens
+
+#ifndef MARS_MOptics
+#include "MOptics.h"
+#endif
+
+class MQuaternion;
+
+class MLens : public MOptics
+{
+private:
+    Double_t fMaxR;
+
+    void InitMaxR();
+
+public:
+    MLens(const char *name=NULL, const char *title=NULL);
+
+    Double_t GetMaxR() const { return fMaxR; }
+    Double_t GetA() const;
+
+    virtual Bool_t CanHit(const MQuaternion &p) const;
+
+    Int_t ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &) const;
+
+    Bool_t IsValid() const { return kTRUE; }
+
+    ClassDef(MLens, 1) // Parameter container storing the description of a lens
+};
+    
+#endif
Index: trunk/Mars/msimreflector/MLens.cc
===================================================================
--- trunk/Mars/msimreflector/MLens.cc	(revision 19592)
+++ 	(revision )
@@ -1,273 +1,0 @@
-/* ======================================================================== *\
-!
-! *
-! * This file is part of CheObs, the Modular Analysis and Reconstruction
-! * Software. It is distributed to you in the hope that it can be a useful
-! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
-! * It is distributed WITHOUT ANY WARRANTY.
-! *
-! * Permission to use, copy, modify and distribute this software and its
-! * documentation for any purpose is hereby granted without fee,
-! * provided that the above copyright notice appears in all copies and
-! * that both that copyright notice and this permission notice appear
-! * in supporting documentation. It is provided "as is" without express
-! * or implied warranty.
-! *
-!
-!
-!   Author(s): Thomas Bretz,  6/2019 <mailto:tbretz@physik.rwth-aachen.de>
-!
-!   Copyright: CheObs Software Development, 2000-2009
-!
-!
-\* ======================================================================== */
-
-//////////////////////////////////////////////////////////////////////////////
-//
-//  MLens
-//
-//////////////////////////////////////////////////////////////////////////////
-#include "MLens.h"
-
-#include <fstream>
-#include <errno.h>
-
-#include "MQuaternion.h"
-
-#include "MLog.h"
-#include "MLogManip.h"
-
-ClassImp(MLens);
-
-using namespace std;
-
-// --------------------------------------------------------------------------
-//
-// Default constructor
-//
-MLens::MLens(const char *name, const char *title)
-{
-    fName  = name  ? name  : "MLens";
-    fTitle = title ? title : "Parameter container storing a collection of several mirrors (reflector)";
-
-    fMaxR = 25;
-}
-
-// --------------------------------------------------------------------------
-//
-// Return the total Area of all mirrors. Note, that it is recalculated
-// with any call.
-//
-Double_t MLens::GetA() const
-{
-    return fMaxR*fMaxR*TMath::Pi();
-}
-
-// --------------------------------------------------------------------------
-//
-// Check with a rough estimate whether a photon can hit the reflector.
-//
-Bool_t MLens::CanHit(const MQuaternion &p) const
-{
-    // p is given in the reflectory coordinate frame. This is meant as a
-    // fast check without lengthy calculations to omit all photons which
-    // cannot hit the reflector at all
-    return p.R2()<fMaxR*fMaxR;
-}
-
-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;
-    }
-};
-
-float 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) * 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)   ( 0 )
-
-   // 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)
-        return 0;
-
-    const double sqrtBac = sqrt(B2 - a*c);
-
-    if (a==0)
-        return 0;
-
-    const double dz[2] =
-    {
-        (sqrtBac-B) / a,
-        (sqrtBac+B) / a
-    };
-
-    // 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;
-}
-
-double CalcRefractiveIndex(const double &lambda)
-{
-    // https://refractiveindex.info/?shelf=organic&book=poly(methyl_methacrylate)&page=Szczurowski
-
-    // n^2-1=\frac{0.99654λ^2}{λ^2-0.00787}+\frac{0.18964λ^2}{λ^2-0.02191}+\frac{0.00411λ^2}{λ^2-3.85727}
-
-    const double l2 = lambda*lambda;
-
-    const double c0 = 0.99654/(1-0.00787/l2);
-    const double c1 = 0.18964/(1-0.02191/l2);
-    const double c2 = 0.00411/(1-3.85727/l2);
-
-    return sqrt(1+c0+c1+c2);
-}
-
-void ApplyRefraction(MQuaternion &u, const TVector3 &n, const double &n1, const double &n2)
-{
-    // u: incoming direction
-    // n: normal vector of surface
-    // n2: refractive index of new(?) medium
-    // n1: refractive index of old(?) medium
-
-    const double r = n2/n1;
-    const double c = n*u.fVectorPart;
-
-    const double R = 1 - r*r*(1-c*c);
-
-    const auto v = r*c + sqrt(R<0 ? 0 : R);
-
-    u = r*u - n*v;
-}
-
-Int_t MLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const
-{
-    const double F = 50;    // [mm] focal length
-    const double R = 25;    // [mm] radius of lens
-    const int    N = 25;    // [cnt] Nuber of grooves within radius
-
-    const double w = R/N;   // [mm] Width of a single groove
-
-    const double n = CalcRefractiveIndex(wavelength==0 ? 450 : wavelength);
-
-    const int nth = TMath::FloorNint(p.R()/w); // Ray will hit the nth groove
-
-    const double r0  =  nth*w;      // Lower radius of nth groove
-    const double r1  = (nth+1)*w;   // Upper radius of nth groove
-    const double cen = (nth+0.5)*w; // Center of nth groove
-
-    // FIXME: Do we have to check the nth-1 and nth+1 groove as well?
-
-    // Angle to the center of the groove
-    const double phi = acos(cen / (F/1.125));
-
-    //        theta                peak_z
-    const Groove g( TMath::Pi()/2 - phi, r0*tan(phi));
-
-    // Calculate the insersection between the ray and the cone
-    float dz = CalcIntersection(p, u, g, -3);
-
-    // No groove was hit
-    if (dz>=0)
-        return -1;
-
-    // Propagate to the hit point
-    p.PropagateZ(u, dz);
-
-    // Check where the ray has hit
-    const double pr = p.R();
-
-    // If the ray has not hit within the right radius.. reject
-    if (pr<=r0 || pr>r1)
-        return -1;
-
-    // Normal vector of the surface at the hit point
-    // FIXME: Surface roughness?
-    TVector3 norm;
-    norm.SetMagThetaPhi(1, g.fTheta+TMath::Pi()/2, p.XYvector().Phi());
-
-    // Apply refraction at lens entrance (change directional vector)
-    // FIXME: Surace roughness
-    ApplyRefraction(u, norm, 1, n);
-
-    // Propagate ray until bottom of lens
-    const double v = u.fRealPart; // c changes in the medium 
-    u.fRealPart = n*v;
-    p.PropagateZ(u, -3);
-    u.fRealPart = v;  // Set back to c
-
-    // New normal surface (bottom of lens)
-    norm.SetMagThetaPhi(1, 0, 0);
-
-    // Apply refraction at lens exit (change directional vector)
-    // FIXME: Surace roughness
-    ApplyRefraction(u, norm, n, 1);
-
-    // 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());
-
-    // This is only correct if the focal distance is calculated from the inner lens surface!
-    p.fVectorPart.SetZ(0);
-
-    return nth; // Keep track of groove index
-}
Index: trunk/Mars/msimreflector/MLens.h
===================================================================
--- trunk/Mars/msimreflector/MLens.h	(revision 19592)
+++ 	(revision )
@@ -1,32 +1,0 @@
-#ifndef MARS_MLens
-#define MARS_MLens
-
-#ifndef MARS_MOptics
-#include "MOptics.h"
-#endif
-
-class MQuaternion;
-
-class MLens : public MOptics
-{
-private:
-    Double_t fMaxR;
-
-    void InitMaxR();
-
-public:
-    MLens(const char *name=NULL, const char *title=NULL);
-
-    Double_t GetMaxR() const { return fMaxR; }
-    Double_t GetA() const;
-
-    virtual Bool_t CanHit(const MQuaternion &p) const;
-
-    Int_t ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &) const;
-
-    Bool_t IsValid() const { return kTRUE; }
-
-    ClassDef(MLens, 1) // Parameter container storing the description of a lens
-};
-    
-#endif