source: trunk/Mars/hawc/fresnellens_psf.C@ 19953

/*****************************************************************

 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_psf()
{
    // ------------------- 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 w = 0.01;   // [cm] Width of a single groove

    // Tim used some unintentionally wrong values in his simulation:
    // const double D = 50.21;   // [cm] Diameter of the refractive surface
    // const double F = 49.833;  // [cm] calculated as D * (F/D) with D=50.21/54.92
    // const double w = 0.05;    // [cm] Width of a single groove

    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 simulated rays

    double psf = 0;  // Do not change! It might have unpredicted side effectes

    double angle = 0; // [cm] angle of incidence of simulates rays

    //TVector3 point_source(0, -F*atan(angle*TMath::DegToRad()), F);  // Point source at x=0, y=0, z=F
    TVector3 point_source; // No point source

    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",     2000*4/3, -R*4/3, R*4/3, 2000, -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;
        if (point_source.Mag()<1e-10)
            dir.SetMagThetaPhi(1, (180-angle)*TMath::DegToRad(), TMath::Pi()/2);
        else
        {
            // Note that this is not a perfect point source as the
            // flux is homogeneous over the surface of the lens
            TVector3 ps = pos - point_source;
            dir.SetMagThetaPhi(1, ps.Theta(), ps.Phi());
        }

        double v = 1./(TMath::C()*100/1e9); // cm/ns

        // ---------------------------------------------------------------

        MQuaternion pp(pos, 0);
        MQuaternion uu(dir, v);

        // ---------------------------------------------------------------

        int ret = lens.ExecuteOptics(pp, uu, lambda);

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