source: trunk/Mars/hawc/fresnellens_image.C@ 19907

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

 fresnellens_image.C - Example how to simulate a fresnel lens

 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_image.C

 or from within root

    [0] .x ../hawc/fresnellens_image.C

******************************************************************/
void fresnellens_image()
{
    Double_t telzd = 0;      // [deg] Telescope orientation Zenith Distance
    Double_t telaz = 0;      // [deg] Telescope orientation Azimuth

    Double_t dzd   = 6;      // [deg] Relative angle of light ray Zd
    Double_t daz   = 0;      // [deg] Relative angle of light ray Az

    Double_t F     = 0.5021; // [m]   Distance of camera to lens (w.r.t. to lens exit)
    //                       // This is different from F in MFresnelLens!

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

    MFresnelLens lens("Reflector");
    //lens.DefineLens(F, D, w, H, lambda);
    //lens.EnableSlopeAbsorption();
    //lens.EnableDraftAbsorption();
    //lens.DisableBottomReflection();
    //lens.DisableMultiEntry();
    //lens.DisableFresnelReflection();
    //lens.SetMinHits(5);
    //lens.SetMaxHits(0);
    lens.ReadTransmission("resmc/hawcseye/transmission-pmma-3mm.txt", 0.3, true);


    MParSpline pde("PhotonDetectionEfficiency");
    MParSpline acc("ConesAngularAcceptance");

    pde.ReadFile("resmc/hawcseye/pde.txt");
    acc.ReadFile("resmc/hawcseye/cones-angular-acceptance.txt");

    MGeomCamFAMOUS geom(F);

    MSimRays simrays;
    simrays.SetNameReflector("Reflector");
    simrays.SetWavelengthRange(300, 890);
    // Use this to simulate a diffent number of photons per bunch
    //simrays.SetNumPhotons(100000);
    // Use this to simulate a point source at distance F form the lens center
    //simrays.SetHeight(F/1000);

    MPointingPos srcpos("Source");
    srcpos.SetLocalPosition(dzd, daz);

    MPointingPos pointpos;
    pointpos.SetLocalPosition(telzd, telaz);

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

    MParList plist;
    MTaskList tlist;
    plist.AddToList(&tlist);

    plist.AddToList(&lens);
    plist.AddToList(&pointpos);
    plist.AddToList(&srcpos);
    plist.AddToList(&geom);
    plist.AddToList(&pde);
    plist.AddToList(&acc);

    tlist.AddToList(&simrays);

    MHPhotonEvent planeG(1, "HPhotonEventGround");
    MHPhotonEvent plane0(2, "HMirrorPlane0");
    MHPhotonEvent plane2(2, "HMirrorPlane2");
    MHPhotonEvent plane3(2, "HMirrorPlane3");
    MHPhotonEvent plane4(2, "HMirrorPlane4");
    MHPhotonEvent planeC(9, "HPhotonEventCamera");
    MHPhotonEvent planeS(9, "HPhotonEventCamera");

    planeG.SetMaxImpact(60);

    MFillH fillG(&planeG, "MPhotonEvent", "FillGround");
    MFillH fill0(&plane0, "MirrorPlane0", "FillReflector");
    MFillH fill2(&plane2, "MirrorPlane2", "FillCandidates");
    MFillH fill3(&plane3, "MirrorPlane3", "FillReflected");
    MFillH fill4(&plane4, "MirrorPlane4", "FillCameraPlane");
    MFillH fillC(&planeC, "MPhotonEvent", "FillCamera");
    MFillH fillS(&planeS, "MPhotonEvent", "FillCamera");

    fillG.SetNameTab("Ground",      "Photon distribution at ground");
    fill0.SetNameTab("LensIn",      "Photon distribution at plane of lens hit surface");
    fill2.SetNameTab("Candidates",  "*Can hit* photon distribution at reflector plane w/ camera shadow");
    fill3.SetNameTab("LensOut",     "Photon distribution at plane output surface");
    fill4.SetNameTab("CameraPlane", "Photon distribution at camera plane");
    fillC.SetNameTab("Camera",      "Photon distribution which hit a cone in the focal plane");
    fillS.SetNameTab("Sensor",      "Photon distribution which hit a sensor in the focal plane");

    MSimAbsorption absapd("SimPhotonDetectionEfficiency");
    MSimAbsorption cones("SimConesAngularAcceptance");
    absapd.SetParName("PhotonDetectionEfficiency");
    cones.SetParName("ConesAngularAcceptance");
    cones.SetUseTheta();
    cones.SetForce();
    absapd.SetForce();
    // We loose another 4% at the cone surface!

    tlist.AddToList(&fillG);
    tlist.AddToList(&absapd);

    MSimReflector simref;
    simref.SetNameReflector("Reflector");
    simref.SetDetectorMargin(-1);
    tlist.AddToList(&simref);

    tlist.AddToList(&fill0);
    tlist.AddToList(&fill2);
    tlist.AddToList(&fill3);
    tlist.AddToList(&fill4);
    tlist.AddToList(&fillC);
    tlist.AddToList(&cones);
    tlist.AddToList(&fillS);

    //tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);

    MStatusDisplay *disp = new MStatusDisplay;

    MEvtLoop loop;
    loop.SetDisplay(disp);
    loop.SetParList(&plist);

    if (!loop.Eventloop(510))
        return;
}