Index: trunk/Mars/hawc/fresnellens_psf.C
===================================================================
--- trunk/Mars/hawc/fresnellens_psf.C	(revision 19645)
+++ trunk/Mars/hawc/fresnellens_psf.C	(revision 19645)
@@ -0,0 +1,133 @@
+/*****************************************************************
+
+ fresnellens_psf.C - Example how to use MFresnelLens directly
+
+ To run the macro from the command line (assuming you are in a directory
+ Mars/build where you have built your Mars environment) you have to do
+
+    root ../hawc/fresnellens_psf.C
+
+ or from within root
+
+    [0] .x ../hawc/fresnellens_psf.C
+
+******************************************************************/
+void fresnellens_psf2()
+{
+    // ------------------- setup lens -----------------------
+
+    const double D = 54.92;  // [cm] Diameter of the refractive surface
+    const double F = 50.21;  // [cm] Nominal focal length of the lens
+    const double H = 0.25;   // [cm] Thickness of lens
+    const double Z = F;      // [cm] camera position w.r.t. lens exit (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 w = 1;      // [cm] Width of a single groove
+    //const double H = 1.5;    // [cm] Thickness of lens
+
+    const double Rcam = 1.5*4.5; // [cm] radius of camera
+
+    // const double lambda = 546; // [nm] Wavelength for lens surface definition
+
+    double psf = 0;  // Do not change! It might have unpredicted side effectes
+
+    double angle = 8;
+
+    MFresnelLens lens;
+    lens.SetPSF(psf);
+    //lens.DefineLens(F, D, w, H, lambda);
+    //lens.EnableSlopeAbsorption();
+    //lens.EnableDraftAbsorption();
+    //lens.DisableBottomReflection();
+    //lens.DisableMultiEntry();
+    //lens.DisableFresnelReflection();
+    lens.ReadTransmission("resmc/hawcseye/transmission-pmma-3mm.txt", 0.3, true);
+
+    // ------------------- setup histograms -----------------
+
+    unsigned int NN = 1000;
+
+    TH2F h_out("OUT", "Photon Distribution (out)",   NN*4/3, -R*4/3, R*4/3,   NN, -R, R);
+    TH2F h_psf("PSF", "Point Spread Function",     1000*4/3, -R*4/3, R*4/3, 1000, -R, R);
+    TH2F h_cam("CAM", "Point Spread Function",     4*37*4/3, -14*4/3, 14*4/3, 4*37, -14, 14);
+
+    TH1F h_rc("RET",  "Return code",                 8, 0.5, 8.5);
+
+    // ------------------ Run simulation ---------------------
+
+    int cnt = 0; // counter for photons hitting the camera
+
+    int Nrays = 100000;
+    for (int i=0; i<Nrays; i++)
+    {
+        // ---------------------------------------------------------------
+        Double_t X, Y;
+        gRandom->Circle(X, Y, R*sqrt(gRandom->Uniform(0, 1)));
+
+        TVector3 pos(X, Y,  0);
+
+        // ---------------------------------------------------------------
+
+        TVector3 dir;
+        dir.SetMagThetaPhi(1, (180-angle)*TMath::DegToRad(), TMath::Pi()/2);
+
+        double v = 1./(TMath::C()*100/1e9); // cm/ns
+
+        // ---------------------------------------------------------------
+
+        MQuaternion pp(pos, 0);
+        MQuaternion uu(dir, v);
+
+        // ---------------------------------------------------------------
+
+        int ret = lens.ExecuteOptics(pp, uu, 546);
+
+        // Skip photons which will not hit the focal plane within r<R
+        if (ret<0)
+        {
+            h_rc.Fill((-ret)%10);  // Keep reason for loss
+            continue;
+        }
+
+        // Keep exit position on bottom surface
+        h_out.Fill(pp.X(), pp.Y());
+
+        // Propagate to th focal plane
+        pp.PropagateZ(uu, Z);
+
+        // Keep position in focal plane
+        h_psf.Fill(pp.X(), pp.Y());
+        h_cam.Fill(pp.X(), pp.Y());
+
+        if (pp.R()<Rcam)
+            cnt++;
+    }
+
+    // ------------------ Display result ---------------------
+
+    TEllipse circ;
+    circ.SetLineWidth(2);
+    circ.SetFillStyle(0);
+
+    TMarker m;
+    m.SetMarkerStyle(kFullDotMedium);
+
+    TCanvas *c = new TCanvas;
+    c->Divide(2,2);
+    c->cd(1);
+    h_cam.DrawCopy("colz");
+    circ.DrawEllipse(0, 0, 4.5*1.5, 4.5*1.5, 0, 360, 0);
+    m.DrawMarker(0, angle*F/TMath::RadToDeg());
+    c->cd(2);
+    h_out.DrawCopy("colz");
+    c->cd(3);
+    h_psf.DrawCopy("colz");
+    circ.DrawEllipse(0, 0, 4.5*1.5, 4.5*1.5, 0, 360, 0);
+    m.DrawMarker(0, angle*F/TMath::RadToDeg());
+    c->cd(4);
+    h_rc.DrawCopy();
+
+    // Photons inside the camera
+    cout << "Efficiency = " << double(cnt)/Nrays << "  (N=" << cnt << ")" << endl;
+}
Index: trunk/Mars/hawc/fresnellens_traceray.C
===================================================================
--- trunk/Mars/hawc/fresnellens_traceray.C	(revision 19645)
+++ trunk/Mars/hawc/fresnellens_traceray.C	(revision 19645)
@@ -0,0 +1,152 @@
+#include <vector>
+#include <iostream>
+
+#include <TH1.h>
+#include <TMarker.h>
+#include <TArrow.h>
+#include <TCanvas.h>
+#include <TRandom.h>
+
+#include "MFresnelLens.h"
+
+/*****************************************************************
+
+ fresnellens_traceray.C - 2D Example of ray-tracing the fresnel lens
+
+ The macro has to be compilesd to work
+
+ To run the macro from the command line (assuming you are in a directory
+ Mars/build where you have built your Mars environment) you have to do
+
+    root ../hawc/fresnellens_traceray.C++
+
+ or from within root
+
+    [0] .x ../hawc/fresnellens_traceray.C++
+
+******************************************************************/
+void fresnellens_traceray()
+{
+    // ------------------- setup lens -----------------------
+
+    const double D = 54.92;  // [cm] Diameter of the refractive surface
+    const double F = 50.21;  // [cm] Nominal focal length of the lens
+    const double Z = F+H/2;  // [cm] camera position w.r.t. lens exit (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] Width of a single groove
+    //const double w = 1;        // [cm] Width of a single groove
+    //const double H = 2.5;    // [cm] Thickness of lens
+
+    const double lambda = 546; // [nm] Wavelength for lens surface definition
+
+    double angle = 0;  // Check ix=0 at angle=-55 and N=10
+
+    double Z0 = 3;
+
+    const double n0 = MFresnelLens::RefractiveIndex(lambda);
+
+    MFresnelLens lens;
+    lens.DefineLens(F, D, w, H, lambda);
+    //lens.SetPSF(psf);
+    //lens.DefineLens(F, D, w, H, lambda);
+    //lens.EnableSlopeAbsorption();
+    //lens.EnableDraftAbsorption();
+    //lens.DisableBottomReflection();
+    //lens.DisableMultiEntry();
+    //lens.DisableFresnelReflection();
+    lens.ReadTransmission("resmc/hawcseye/transmission-pmma-3mm.txt", 0.3, true);
+
+    // ------------------ Setup display ---------------------
+
+    TH1C h("", "", 10000, -R*1.1, R*1.1);
+    h.SetStats(kFALSE);
+    h.SetMaximum( 3);
+    h.SetMinimum(-1.1*F);
+    h.DrawCopy();
+
+    TLine line;
+    line.SetLineWidth(2);
+
+    for (unsigned int i=0; i<lens.GetNumGrooves(); i++)
+    {
+        const MFresnelLens::Groove &g = lens[i];
+
+        line.SetLineColor(kBlack);
+        line.DrawLine( w*i, 0,  g.r, g.slope.h);
+        line.DrawLine(-w*i, 0, -g.r, g.slope.h);;
+
+        line.SetLineColor(kRed);
+        line.DrawLine( g.r, g.slope.h,  w*i+w, 0);
+        line.DrawLine(-g.r, g.slope.h, -w*i-w, 0);
+    }
+
+    line.SetLineColor(kBlack);
+    line.SetLineWidth(1);
+    line.DrawLine(-30, -H,   30, -H);
+    line.DrawLine(-30, -Z, 30, -Z);
+    line.DrawLine( R,  -Z,  R,  0);
+    line.DrawLine(-R,  -Z, -R,  0);
+    line.SetLineWidth(3);
+    line.DrawLine(-4.5*1.5, -Z, 4.5*1.5, -Z);
+    line.SetLineWidth(1);
+
+    // ------------------------------------------------------
+
+    TArrow arrow;
+    TMarker marker;
+
+    marker.SetMarkerColor(kBlue);
+    marker.SetMarkerStyle(kFullDotMedium);
+
+    arrow.SetArrowSize(0.01);
+    arrow.SetLineColor(kBlue);
+
+    double vc = 1./(TMath::C()*100/1e9); // cm/ns
+
+    for (int i=0; i<100; i++)
+    {
+        double X = gRandom->Uniform(-R, R);
+
+        TVector3 pos(X, 0,  0);
+        TVector3 dir(0, 0, -1);
+
+        dir.SetMagThetaPhi(1, (180-angle)*TMath::DegToRad(), 0);
+
+        MQuaternion p(pos, 0);
+        MQuaternion u(dir, vc);
+
+        // Photon must be at the lens surface before it can enter the lens
+        p.PropagateZ(u, Z0);
+
+        vector<MQuaternion> vec;
+        vec.push_back(p);
+
+        const int cnt = lens.TraceRay(vec, p, u, lambda, true);
+
+        // Particle sucdessfully reached the focal plane
+        if (cnt>=0)
+        {
+            p.PropagateZ(u, -Z);
+            vec.push_back(p);
+
+            p.PropagateZ(u, u.Z()<0 ? -2.5*F : Z0);
+            vec.push_back(p);
+        }
+
+        // Particle is upgoing and has hit no surface
+        if (cnt==-1291 || cnt==-1391)
+        {
+            p.PropagateZ(u, 3);
+            vec.push_back(p);
+        }
+
+        arrow.SetLineColor(kBlue);
+        for (unsigned int j=0; j<vec.size()-1; j++)
+        {
+            arrow.DrawArrow(vec[j].X(), vec[j].Z(), vec[j+1].X(), vec[j+1].Z());
+            marker.DrawMarker(vec[j+1].X(), vec[j+1].Z());
+        }
+    }
+}