Modules/FdSimulation/TelescopeSimulatorKG2/Camera.cc

Go to the documentation of this file.

//#define DEBUG

#include "Camera.h"
#include "RTFunctions.h"
#include "RayTracer.h"

#include <utl/MathConstants.h>
#include <utl/AugerUnits.h>
#include <utl/ErrorLogger.h>
#include <utl/Photon.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/AugerException.h>
#include <utl/RandomEngine.h>
#include <utl/Math.h>

#include <CLHEP/Random/Randomize.h>

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Camera.h>
#include <fdet/Pixel.h>
#include <fdet/Channel.h>
#include <fdet/Mirror.h> // for check of photons inside camera bay

#include <TPolyLine3D.h>
#include <TLine.h>
#include <TMarker.h>
#include <TLatex.h>
#include <TPolyMarker3D.h>
#include <TCanvas.h>
#include <TView.h>
#include <TPad.h>
#include <TObjArray.h>

#include <sstream>
#include <iostream>
#include <iomanip>
#include <vector>
#include <utility>

using namespace std;
using namespace TelescopeSimulatorKG2;
using namespace utl;

using CLHEP::RandFlat;


Camera::Camera(RandomEngine& rndm,
               const fdet::Telescope& tel,
               const bool plotPhotonTracksAtMercedes,
               const bool simulateMercedesStars,
               const bool simulateHaloEffects,
               const bool simulateGhostEffects,
               const double pmt_n) :
  fRandom(rndm),
  fTel(tel),
  fPlotPhotonTracksAtMercedes(plotPhotonTracksAtMercedes),
  fSimulateMercedesStars(simulateMercedesStars),
  fSimulateHaloEffects(simulateHaloEffects),
  fSimulateGhostEffects(simulateGhostEffects),
  fPMT_n(pmt_n),
  fTelCS(tel.GetTelescopeCoordinateSystem()),
  fOrigin(0, 0, 0, fTelCS),
  fRFocal(tel.GetCamera().GetRadiusFocal()),
  fMercedesEff(tel.GetCamera().GetMercedesEfficiency()),
  fSigma(tel.GetCamera().GetSigmaNormal()),
  fRCurv2(Sqr(tel.GetMirror().GetRadiusOfCurvature()))  // need this only to check for photons to stay in mirror bay
{
  // Camera SetMercedesParameter
  // Find the position of the vertex in the reference frame of the
  // hexagon. See picture.
  /*
          ^      /\
         n \    /  \
            \  /    \
             \/      \     height
             /        \
        ^   /          \
       z|  /            \
        |  --------------
        |       2xwidth
    <----
     t
  */
  // Set mercedes parameters
  const fdet::Camera& camera = tel.GetCamera();
  fHMerc = camera.GetMercedesHeight();  // height of mercedes ~ 17.91 mm 

  const double ddmerc = camera.GetMercedesBase();  // ~ 4.8 mm = half base length (full 9.2mm)
  const double dTheta = 1.5*deg; // For numbers see e.g. GAP1999_027
  const double dPhi = dTheta * cos(30*deg); // ~1.3*deg
  fR2 = fRFocal * sin(2 * dPhi / 3); // pixel side length ~ 26.34*mm;
  fD2 = fRFocal * sin(dTheta); // horizontal full pixel width ~ 45.6*mm
  fD1 = fD2 - 2 * ddmerc;    // horizontal pixel width (inside mercedes)
  fR1 = fD1 / utl::kSqrt3;  // hexagon side length (inside mercedes) [Srqt(3) = 2 * cos(30*deg)]
  fNz = 0.5 * (fD2 - fD1);  // = ddmerc; mercedes width
  fNt = fHMerc;             // mercedes height
  const double aux = sqrt(fNz * fNz + fNt * fNt);
  fNz /= aux;               // (Nz, Nt) : normal on mercedes CHECK !! Seems to be wrong
  fNt /= aux;

  //CalculateCameraSupportCoordinateSystem
  const double theta = 16*deg;
  const double shift = fRFocal - 50*cm;
  //const double shift = fRFocal - 30 * cm;
  const Vector vShift(0, 0, -1, fTelCS);
  const CoordinateSystemPtr& auxCS = CoordinateSystem::Translation(vShift * shift, fTelCS);
  fCameraSupportCoordinateSystem = CoordinateSystem::RotationY(theta, auxCS);
}


RTNext
Camera::Trace(const Photon& photonIn,
              Photon& photonOut,
              int& nreflections,
              int& nbackscattered,
              int& col,
              int& row)
{
  // simulate camera, including/excluding the mercedes

  const bool backwards = photonIn.GetDirection().GetZ(fTelCS) < 0 ; // towards the mirror 

  utl::Photon photonCamera;
  if (!TraceSurf(photonIn, photonCamera)) {
    // missed camera sphere
    photonOut = photonIn;
    return fSimulateGhostEffects ? backwards ? eMirror : eCameraShadow : eLost;
  }

  // theoretical incidence direction of the photon - use just to find the pixelId 
  const Vector incidentPhotonDir = Normalized(photonCamera.GetPosition() - fOrigin);
  if (!FindPixelId(incidentPhotonDir, row, col)) {
    // camera missed ...
    photonOut = photonCamera;
    return fSimulateGhostEffects ? backwards ? eMirror : eCameraShadow : eLost;
  }

  const CoordinateSystemPtr& pxlCS = fTel.GetPixel(row, col).GetPixelCoordinateSystem();

  if (fSimulateMercedesStars) {
    // simulate mercedes
    const RTResult mercStatus = TraceMerc(pxlCS, photonCamera, photonOut, nreflections);
    if (mercStatus == eInvalidMercedes)
      return eDrawAndExit;

    if (fSimulateHaloEffects &&
        (mercStatus == eBackscattered ||
         mercStatus == eInvalidHitMercedes)) {
      ++nbackscattered;
      return eMirror;
    }

    if (mercStatus == eOK)
      return eDetected;

  } else {  // no mercedes

    // measurements in 2010 showed that the PMT reflectivity is ~20%
    // the PMT refractive index can be set in the TelescopeSimulatorKG2.xml.in, default is 2.58 -> R=19.5%
    const double r_index_12 = fAir_n / fPMT_n;
    RTFunctions::PhotonList photonsHittingPhotocathode;
    const Vector normalPMT(0, 0, -1, pxlCS);

    RTFunctions::Refraction(r_index_12,
                            photonCamera,
                            normalPMT,
                            photonsHittingPhotocathode);
    const utl::Photon& photonRefracted = photonsHittingPhotocathode[0];
    // photonsHittingPhotocathode[0] ist refracted
    // photonsHittingPhotocathode[1] ist reflected

    if (fSimulateHaloEffects) {

      const double tmpRandNr = RandFlat::shoot(&fRandom.GetEngine(), 0, 1);

      // photonRefracted.GetWeight() = photonPMT.GetWeight() * T
      // => T = photonRefracted.GetWeight()/photonPMT.GetWeight()

      const double weight = photonCamera.GetWeight();
      if (tmpRandNr > photonRefracted.GetWeight() / weight) {
        // photon bunch is considered as reflected
        photonOut = photonsHittingPhotocathode[1];
        photonOut.SetWeight(weight);
        // debug
        //cout << "reflected at camera (no mercedes)" << endl;
        return eMirror;
      } else {
        // photon bunch is refracted
        photonOut = photonRefracted;
        photonOut.SetWeight(weight);
        return eDetected;
      }
    }

    // photonOut is not correct (should be photonRefracted) but was
    // implemented like this in the previous TelescopeSimulatorKG
    // -> remove some day
    photonOut = photonRefracted;  // the reflected component is killed at this point, the weight of reflected photons is lost
    return eDetected;  // absorbed by photo-cathode
  } // end: no mercedes

  return eLost;
}


bool
Camera::FindPixelId(const Vector& incidentDir, int& row, int& col) const
{
  const double fNx0 = 10;
  const double fNy0 = 35 / 3; //=11.66666666666667;  // optical axis
  const double dOmega = fTel.GetCamera().GetEta();
  const double dPhi = 0.5 * kSqrt3 * dOmega;
  const int nCols = fTel.GetLastColumn() - fTel.GetFirstColumn() + 1;  // 20
  const int nRows = fTel.GetLastRow() - fTel.GetFirstRow() + 1;  // 22

  const double x = incidentDir.GetX(fTelCS);
  const double y = incidentDir.GetY(fTelCS);
  const double z = incidentDir.GetZ(fTelCS);

  const double omegap = asin(y);
  const double phip = atan2(-x, -z);

  // find nx and ny from angles omegap, phip.
  // -------------------------------------------------------------------
  // ------ algorit. developed by Ronald Shellard, in the FDSim --------
  // ------------------ (FORTRAN version) ------------------------------
  // -------------------------------------------------------------------
  double nxl = 2 * (omegap / dOmega) + 2 * fNx0;
  double nyl = fNy0 - phip / dPhi;

  if (nxl < 0 || nxl > 2 * nCols + 1 ||
      nyl < 0 || nyl > nRows + 1)
    return false;  // the ray of light does not reach any pmt

  col = int(nxl);
  row = int(nyl);
  nxl -= col;
  nyl -= row;

  const double f1 = 3 * nyl - nxl - 1;
  const double f2 = 3 * nyl + nxl - 2;
  double f = 0;

  if (col % 2) {  // col is odd
    col = (col + 1) / 2;
    f = (row % 2) ? f1 : f2;
  } else {
    if (!(row % 2)) {
      f = f1;
      col = col / 2 + (f > 0 ? 0 : 1);
    } else {
      f = f2;
      col = col / 2 + (f > 0 ? 1 : 0);
    }
  }

  if (f >= 0)
    ++row;

  return col >= 1 && row >= 1 && col <= nCols && row <= nRows;
}


RTNext
Camera::TraceShadow(const utl::Photon& photonIn,
                    utl::Photon& photonOut,
                    const bool simulateCameraSupport)
{
  // simulate shadow effect of the camera

  const bool backwards = photonIn.GetDirection().GetZ(fTelCS) > 0 ;  // towards the aperture !!

  RTFunctions::IntersectionList intersections;

  const bool debug = false;

  // Warning: the values of SIN_?? may need some changes to simulate
  // correctly the camera (shadow and trace_surf).
  // side, top and down limits
/*
  const double fSinAx = 0.2491262;
  const double fSinBx = 0.25646528973;
  const double fSinAy = 0.27966724;
  const double fSinBy = 0.26704475;
  const double rCam = fRFocal;
*/

  bool absorbedByCameraBody = false;
  const double rCamIn = 1.641*m;  // see GAP_1999_027
  Point pp;
  if (RTFunctions::Sphere(fOrigin, rCamIn, photonIn, intersections)) {
    if (intersections.size() != 2)
      return backwards ? eLens : eMirror;  
    // select the intersection point with NEGATIVE z-coordinate -> location of camera !
    const bool firstIntersection = intersections.front().first.GetPosition().GetZ(fTelCS) < 0;
/*
    const Vector normal(firstIntersection ? intersections.front().second : intersections.back().second);
    const double x = normal.GetX(fTelCS);
    const double y = normal.GetY(fTelCS);

    if (-fSinAx < x && x < fSinBx &&
        -fSinAy < y && y < fSinBy) {
      photonOut = firstIntersection ? intersections.front().first : intersections.back().first;
      return eLost;  // the light hit the PMT camera back
    }
*/
    Photon p(firstIntersection ? intersections.front().first : intersections.back().first);
    pp = p.GetPosition();
    const double x = abs(pp.GetX(fTelCS));  // symetric
    const double y = abs(pp.GetY(fTelCS));  // symetric
    if (x <= (0.86*m / 2) && y <= (0.93*m / 2))
      absorbedByCameraBody = true;
  }

/*
  const double rCamOut = 1.701*m;  // see GAP_1999_027
  if (!absorbedByCameraBody && RTFunctions::Sphere(fOrigin, rCamOut, photonIn, intersections)) {
    if (intersections.size() != 2)
      return backwards ? eLens : eMirror;  
    const bool firstIntersection = intersections.front().first.GetPosition().GetZ(fTelCS) < 0;
    Photon p(firstIntersection ? intersections.front().first : intersections.back().first);
    pp = p.GetPosition();
    const double x = abs(pp.GetX(fTelCS));  // symetric
    const double y = abs(pp.GetY(fTelCS));  // symetric
    if (x <= (0.86*m / 2) && y <= (0.93*m / 2))
      absorbedByCameraBody = true;
  }
*/

  if (absorbedByCameraBody) {
    if (debug)
      cout << pp.GetX(fTelCS)/mm << ',' << pp.GetY(fTelCS)/mm << ',' << pp.GetZ(fTelCS)/mm << '\n';
    return eLost;  // the light hit the camera back
  }


  if (simulateCameraSupport) {
    const Point pCamera(0,0,0, fCameraSupportCoordinateSystem);
    const Vector nSupport(0,0,1, fCameraSupportCoordinateSystem);
    bool absorbedByCamSupport = false;
    if (RTFunctions::Plane(pCamera, nSupport, photonIn, photonOut)) {
      const Point& pPlane = photonOut.GetPosition();
      const double x = pPlane.GetX(fCameraSupportCoordinateSystem);
      const double y = abs(pPlane.GetY(fCameraSupportCoordinateSystem));  // symetric
      // check if photons on camera-support-plane are still within bay (within mirror sphere)
      const double rPhoton2 = (pPlane - fOrigin).GetMag2();
      if (rPhoton2 > fRCurv2) 
        return eLost;
      if (x < 0)  // camera support is only below the camera !
        return backwards ? eLens : eMirror; 
      
      const double ySupportBegin = 35*cm;
      const double ySupportEnd = ySupportBegin + 5*cm;
      const double zSupportDepth = 10*cm;
      if (y >= ySupportBegin && y <= ySupportEnd) 
        absorbedByCamSupport = true;

      /*
      Each leg has the size of 100x50mm ( C-shaped )
      Longer side is parallel to telescope axis! (See GAP-1999-027)
      */
      const Vector nLongerSideSupport(0, 1, 0, fCameraSupportCoordinateSystem);
      Photon dummyPhoton = photonOut;
      // right bar left
      const Point pLongerSideSupport1(0, ySupportBegin, 0, fCameraSupportCoordinateSystem);
      if (!absorbedByCamSupport && RTFunctions::Plane(pLongerSideSupport1, nLongerSideSupport, photonIn, dummyPhoton)) {
        const double z = dummyPhoton.GetPosition().GetZ(fCameraSupportCoordinateSystem);
        if (z <= 0*cm && z >= -zSupportDepth)
          absorbedByCamSupport = true;
      }
      // right bar right
      const Point pLongerSideSupport2(0, ySupportEnd, 0, fCameraSupportCoordinateSystem);
      if (!absorbedByCamSupport && RTFunctions::Plane(pLongerSideSupport2, nLongerSideSupport, photonIn, dummyPhoton)) {
        const double z = dummyPhoton.GetPosition().GetZ(fCameraSupportCoordinateSystem);
        if (z <= 0*cm && z >= -zSupportDepth)
          absorbedByCamSupport = true;
      }
      // left bar right
      const Point pLongerSideSupport3(0, -ySupportBegin, 0, fCameraSupportCoordinateSystem);
      if (!absorbedByCamSupport && RTFunctions::Plane(pLongerSideSupport3, nLongerSideSupport, photonIn, dummyPhoton)) {
        const double z = dummyPhoton.GetPosition().GetZ(fCameraSupportCoordinateSystem);
        if (z <= 0*cm && z >= -zSupportDepth)
          absorbedByCamSupport = true;
      }
      // left bar left
      const Point pLongerSideSupport4(0, -ySupportEnd, 0, fCameraSupportCoordinateSystem);
      if (!absorbedByCamSupport && RTFunctions::Plane(pLongerSideSupport4, nLongerSideSupport, photonIn, dummyPhoton)) {
        const double z = dummyPhoton.GetPosition().GetZ(fCameraSupportCoordinateSystem);
        if (z <= 0*cm && z >= -zSupportDepth)
          absorbedByCamSupport = true;
      }
      if (absorbedByCamSupport) {
        if (debug) {
          const Point p = dummyPhoton.GetPosition();
          cout << p.GetX(fTelCS)/mm << ',' << p.GetY(fTelCS)/mm << ',' << p.GetZ(fTelCS)/mm << '\n';
        }
        return eLost;
      }
    }
  } else { // do not simulate the camera support
    const bool firstIntersection = intersections.front().first.GetPosition().GetZ(fTelCS) < 0;
    photonOut = firstIntersection ? intersections.front().first : intersections.back().first;
  }
  return backwards ? eLens : eMirror;
}


// find the point that the light hit the focal surface
bool
Camera::TraceSurf(const utl::Photon& photonIn, utl::Photon& photonOut)
{
  RTFunctions::IntersectionList intersections;
  if (RTFunctions::Sphere(fOrigin, fRFocal, photonIn, intersections)) {
    photonOut = intersections.front().first;
    return true;  // hit the camera sphere
  }
  return false;  // miss the camera sphere
}


// simulate Mercedes/PMT system
RTResult
Camera::TraceMerc(const CoordinateSystemPtr& pxlCS,
                  const utl::Photon& photonIn, utl::Photon& photonOut,
                  int& nreflections)
{
  //const Vector& dir = photonIn.GetDirection();

  //const Point pp = photonIn.GetPosition();
  //cout << pp.GetX(fTelCS)/mm << ',' << pp.GetY(fTelCS)/mm << ',' << pp.GetZ(fTelCS)/mm << '\n';

  utl::Photon photon = photonIn;  // starting point

#ifdef DEBUG
  cout << "TraceMerc: photon-in"
          " x=" <<  setprecision(3) << setw(6) << photon.GetPosition().GetX(pxlCS)
       << " y=" <<  setprecision(3) << setw(6) << photon.GetPosition().GetY(pxlCS)
       << " z=" <<  setprecision(3) << setw(6) << photon.GetPosition().GetZ(pxlCS)
       << " nx=" <<  setprecision(3) << setw(6) << photon.GetDirection().GetX(pxlCS)
       << " ny=" <<  setprecision(3) << setw(6) << photon.GetDirection().GetY(pxlCS)
       << " nz=" <<  setprecision(3) << setw(6) << photon.GetDirection().GetZ(pxlCS)
       << endl;
#endif

  // simulate reflections at the mercedes
  utl::Photon point_next;  // the next point of ray-tracing
  int iPlanePrev = 7;  // coming from focal plane !

  do {  // loop until the light hits photo cathode

    bool firstPlane = true;
    double closestDist = 0;
    int closestPlane = 0;
    Vector nClosestPlane;

    const double eps = 0*mm;

    // find the closest plane that light hits (cathode, six merce-walls, focal plane)
    for (int iPlane = 0; iPlane <= 7; ++iPlane) {

      if (iPlane == iPlanePrev)
        continue;

      // 0 to photo cathode surface and 1,2,3,4,5 or 6 to mercedes, 7 focal plane (pixel exit)

#warning cache planes for better performance
      Point pPlane;
      Vector nPlane;
      ConstructPlane(iPlane, nPlane, pPlane, pxlCS);

      Photon testPhoton;
      const double dist = RTFunctions::Plane(pPlane, nPlane, photon, testPhoton);

      if (dist < -eps)
        continue;

      if (firstPlane) {
        firstPlane = false;
        closestDist = dist;
        closestPlane = iPlane;
        nClosestPlane = nPlane;
        point_next = testPhoton;

      } else if (dist < closestDist) {

        closestDist = dist;
        closestPlane = iPlane;
        nClosestPlane = nPlane;
        point_next = testPhoton;

      }
    } // plane loop

    // debug
    //cout << "closestPlane = " << closestPlane << endl;

#ifdef DEBUG
    cout << "TraceMerc: nmerc=" << setw(3) << nreflections
         << " x=" <<  setprecision(3) << setw(6) << point_next.GetPosition().GetX(pxlCS)
         << " y=" <<  setprecision(3) << setw(6) << point_next.GetPosition().GetY(pxlCS)
         << " z=" <<  setprecision(3) << setw(6) << point_next.GetPosition().GetZ(pxlCS)
         << " nx=" <<  setprecision(3) << setw(6) << point_next.GetDirection().GetX(pxlCS)
         << " ny=" <<  setprecision(3) << setw(6) << point_next.GetDirection().GetY(pxlCS)
         << " nz=" <<  setprecision(3) << setw(6) << point_next.GetDirection().GetZ(pxlCS)
         << " " << iPlanePrev << " -> " << closestPlane
         << endl;
#endif

    iPlanePrev = closestPlane;

    if (fPlotPhotonTracksAtMercedes && fPhotonTrack)
      fPhotonTrack->AddPoint(point_next.GetPosition());

    if (firstPlane) {
      ostringstream err;
      err << " RT::Camera>  NO MERCEDES OR PMT HIT! PHOTON invalid! w=" << photonIn.GetWeight()
          << " closestDist=" << closestDist/mm << " mm";
      ERROR(err);
      photonOut = photonIn;
      photonOut.SetPosition(photonOut.GetPosition() + photonOut.GetDirection()*0.5*m);
      return eInvalidHitMercedes;
    }

    if (closestPlane == 0) {  // hit PMT

      // measurements in 2010 showed that the PMT reflectivity is ~20%
      // the PMT refractive index can be set in the TelescopeSimulatorKG2.xml.in, default is 2.58 -> R=19.5%
      const double r_index_12 = fAir_n / fPMT_n;
      RTFunctions::PhotonList photonsHittingPhotocathode;
      RTFunctions::Refraction(r_index_12, point_next, nClosestPlane, photonsHittingPhotocathode);
      const utl::Photon& photonRefracted = photonsHittingPhotocathode[0];
      // photonsHittingPhotocathode[0] ist refracted
      // photonsHittingPhotocathode[1] ist reflected

      if (false) { // photons get absorbed by photo-cathod
        // for estimating the camera surface backreflection
        photonOut = photonRefracted;
        return eOK;  
      }

      if (fSimulateHaloEffects) {

        const double tmpRandNr = RandFlat::shoot(&fRandom.GetEngine(), 0, 1);

        // photonRefracted.GetWeight() = photonPMT.GetWeight() * T
        // => T = photonRefracted.GetWeight()/photonPMT.GetWeight()

        const double weight = point_next.GetWeight();
        if (tmpRandNr > photonRefracted.GetWeight() / weight) {
          // photon bunch is considered as reflected
          photon = photonsHittingPhotocathode[1];
          photon.SetWeight(weight);  // backscattered from cathode
          // diffusive component
          //double amountOfStrayLight = 0.05; // 5% of the light will go in "lambertian directions"
          double amountOfStrayLight = 0;  // 0% of the light will go in "lambertian directions"
          if (amountOfStrayLight > 0) {
            const double maxTheta = 25 * degree;
            const Vector lambertdiffusion =
              RTFunctions::LambertDiffusion(fRandom, photon.GetDirection(), nClosestPlane, amountOfStrayLight, maxTheta);
            photon.SetDirection(lambertdiffusion);
          }
        } else {
          // photon bunch is refracted
          photonOut = photonRefracted;
          photonOut.SetWeight(weight);
          return eOK;
        }

      } else {  // no halo:

        photonOut = photonRefracted;  // the reflected component is killed at this point, the weight of reflected photons is lost
        return eOK;  // absorbed by photo-cathode
      }

    } else if (closestPlane == 7) {  // hit focal surface (upper range of mercedes, leaving pixel towards mirror)

      // only when we hit the surface from the INSIDE
      if (point_next.GetDirection() * nClosestPlane <= 0) {
        photonOut = point_next;
        if (photonOut.GetDirection().GetZ(fTelCS) > 0) {
          ostringstream err;
          err << "Photon backscattered from focal surface towards aperture! ";
          INFO(err);
        }
        return eBackscattered;
      }

      photon = point_next;

    } else {  // -> We hit a mercedes

      RTFunctions::Reflection(point_next,
                              nClosestPlane,  // normal vector of plane
                              photon);  // reflected photon
      photon.SetWeight(point_next.GetWeight() * fMercedesEff);
      // diffusive component
      //double amountOfStrayLight = 0.05; // 5% of the light will go in "lambertian directions"
      double amountOfStrayLight = 0; // 0% of the light will go in "lambertian directions"
      if (amountOfStrayLight > 0 ) {
        const double maxTheta = 25*degree;
        const Vector lambertdiffusion =
          RTFunctions::LambertDiffusion(fRandom, photon.GetDirection(), nClosestPlane, amountOfStrayLight, maxTheta);
        //const Vector lambertdiffusion = RTFunctions::LambertDiffusion( fRandom, photon.GetDirection(), photon.GetDirection(), amountOfStrayLight, maxTheta);
        //cout << "carefull!!!!!, amountOfStrayLight = 45%, maxTheta = 5degr, distribution peaked around outgoing photon direction (instead of surface normal)" << endl;
        photon.SetDirection(lambertdiffusion);
      }

      ++nreflections;  // Do next reflection
      /* // mercedes black
      #warning mercedes is black -> photons should be lost, here weight = 0
      photon.SetWeight(0.);
      */

    }

  } while (true);  // WHILE we hit the PMT
}


/*
              /\
             /  \
            /    \
          R/      \R
          /        \
         /  2    1  \
        |            |
        |            |
       R|  3   .   6 |R
        |            |
        |            |
         \   4   5  /
          \        /
          R\      /R
  Y^        \    /
   |         \  /
   |          \/
   |____>
        X

   R: length of vertex
 */

Camera::HexagonDirection
Camera::Hexagon(const double x, const double y, const double r)
{
  // verify if position x,y is inside or outside of a hexagon,
  // which the distance from the center (0,0) to the vertex
  // is smaller than r (see picture)
  // ***************
  // 0 -> inside the hexagon
  // from 1 to 6 -> outside the hexagon

  const double xl = std::fabs(x);
  const double yl = std::fabs(y);
  const double xm = 0.5 * utl::kSqrt3 * r;  // cos(30)*r
  const double tmp = 1 / utl::kSqrt3;  // tan(30)

  // inside hexagon check
  if (xl <= xm && yl <= r - tmp * xl)
    return eInside;

  if (yl < tmp * xl)  // left and right slice
    return (x > 0) ? eRight : eLeft;
  else  // upper right half
    return (y > 0) ? (x > 0 ? eTopRight:eTopLeft) : (x > 0 ? eBottomRight:eBottomLeft);
}


void
Camera::ConstructPlane(const int id, Vector& n, Point& p, const CoordinateSystemPtr& pxlCS)
{
  // define the plane surfaces to simulate the mercedes
  // --------------------------
  // id=0 -> photo-cathode surface
  //    7 -> focal surface
  //    1,2,3,4,5 and 6 -> mercedes surfaces (see picture)
  // ++++++++++++++++++++++++++
  // n -> normal vector of the surface
  // o -> a point of the plane surface
  // *************************
  // 0 -> there is no plane identified by id
  // 1 -> ok

  if (id == 0 || id == 7) {
    p = Point(0, 0, id == 0 ? 0 : fHMerc, pxlCS);
    n = Vector(0, 0, id == 0 ? 1 : -1 , pxlCS);
    return;
  }

  if (id > 0 && id < 7) {
    const double dTheta = 60 * deg;
    const double theta = dTheta * id;
    const double sinTheta = sin(theta);  // y
    const double cosTheta = cos(theta);  // x
    n = Vector(-fNt * cosTheta, -fNt * sinTheta, fNz, pxlCS);
    p = Point(0.5 * cosTheta * fD1, 0.5 * sinTheta * fD1, 0, pxlCS);
    return;
  }

  ostringstream err;
  err << "This can never happen! Fix TelescopeSimulator. There is no such plane as defined by id=" << id;
  ERROR(err);
  throw utl::AugerException(err.str());
}


TObjArray*
Camera::Draw(const bool do3D)
{
  TObjArray* const objs = new TObjArray();

  double xMin = 0;  // <-- TODO: replace with utl::MinMax<double>
  double xMax = 0;
  double yMin = 0;
  double yMax = 0;
  double zMin = 0;
  double zMax = 0;
  bool first = true;

  const int nCols = fTel.GetLastColumn() - fTel.GetFirstColumn() + 1;  // 20
  const int nRows = fTel.GetLastRow() - fTel.GetFirstRow() + 1;  // 22

  // int color = 0;
  for (int col = 1; col <= nCols; ++col) {

    for (int row = 1; row <= nRows; ++row) {

      const int pixelId = (col - 1) * nRows + row;
      const CoordinateSystemPtr& pxlCS = fTel.GetPixel(row, col).GetPixelCoordinateSystem();

      const double r_in = fR1;
      const double r_out = fR2;

      for (int i = 0; i < 6; ++i) {

        const double theta1 = 30*deg + 60*deg * i;
        const double theta2 = theta1 + 60*deg;

        const double s1 = sin(theta1);
        const double c1 = cos(theta1);
        const Point p1_in(c1 * r_in, s1 * r_in, 0, pxlCS);

        const double s2 = sin(theta2);
        const double c2 = cos(theta2);
        const Point p2_in(c2 * r_in, s2 * r_in, 0, pxlCS);

        const Point p1_out(c1 * r_out, s1 * r_out, fHMerc, pxlCS);
        const Point p2_out(c2 * r_out, s2 * r_out, fHMerc, pxlCS);

        // pixel position in telCS
        const Point zeroPx(0,0,0, pxlCS);
        const double pxl_x = zeroPx.GetX(fTelCS);
        const double pxl_y = zeroPx.GetY(fTelCS);
        const double pxl_z = zeroPx.GetZ(fTelCS);

        if (first) {
          first = false;
          xMin = pxl_x;
          xMax = pxl_x;
          yMin = pxl_y;
          yMax = pxl_y;
          zMin = pxl_z;
          zMax = pxl_z;
        } else {
          if (pxl_x > xMax)
            xMax = pxl_x;
          if (pxl_y > yMax)
            yMax = pxl_y;
          if (pxl_z > zMax)
            zMax = pxl_z;
          if (pxl_x < xMin)
            xMin = pxl_x;
          if (pxl_y < yMin)
            yMin = pxl_y;
          if (pxl_z < zMin)
            zMin = pxl_z;
        }

        if (do3D) {

          /*TPolyMarker3D *p_pos = new TPolyMarker3D(1);
          p_pos->SetPoint(0, pxl_x, pxl_y, pxl_z);
          p_pos->SetMarkerColor(color++);
          p_pos->SetMarkerStyle(20);
          p_pos->SetMarkerSize(.3);
          p_pos->Draw();
          objs->AddLast(p_pos);*/

        } else {
          //double w, h;
          ostringstream pxl_txt;
          pxl_txt << pixelId;
          TLatex* const p_pos = new TLatex(pxl_x, pxl_y, pxl_txt.str().c_str());
          p_pos->SetX(pxl_x - p_pos->GetXsize() / 6.);
          p_pos->SetY(pxl_y - p_pos->GetYsize() / 6.);
          p_pos->SetTextColor(1);
          p_pos->SetTextSize(0.01);
          p_pos->SetTextAngle(0);
          p_pos->Draw();
          objs->AddLast(p_pos);
        }

        if (do3D) { // DO THE 3D PLOT

          TPolyLine3D* const l_in = new TPolyLine3D(1);
          l_in->SetLineColor(2);
          l_in->SetPoint(0, p1_in.GetX(fTelCS), p1_in.GetY(fTelCS), p1_in.GetZ(fTelCS));
          l_in->SetPoint(1, p2_in.GetX(fTelCS), p2_in.GetY(fTelCS), p2_in.GetZ(fTelCS));
          l_in->Draw();
          objs->AddLast(l_in);

          TPolyLine3D* const l_out = new TPolyLine3D(1);
          l_out->SetLineColor(1);
          l_out->SetLineStyle(2);

          l_out->SetPoint(0, p1_out.GetX(fTelCS), p1_out.GetY(fTelCS), p1_out.GetZ(fTelCS));
          l_out->SetPoint(1, p2_out.GetX(fTelCS), p2_out.GetY(fTelCS), p2_out.GetZ(fTelCS));
          l_out->Draw();
          objs->AddLast(l_out);

          TPolyLine3D* const l_corner = new TPolyLine3D(1);
          l_corner->SetLineColor(1);
          l_corner->SetLineStyle(2);

          l_corner->SetPoint(0, p1_in.GetX(fTelCS), p1_in.GetY(fTelCS), p1_in.GetZ(fTelCS));
          l_corner->SetPoint(1, p1_out.GetX(fTelCS), p1_out.GetY(fTelCS), p1_out.GetZ(fTelCS));
          l_corner->Draw();
          objs->AddLast(l_corner);

        } else { // DO THE 2D PLOT

          TLine* const l_in = new TLine(p1_in.GetX(fTelCS), p1_in.GetY(fTelCS),
                                        p2_in.GetX(fTelCS), p2_in.GetY(fTelCS));
          l_in->SetLineColor(2);

          TLine* const l_out = new TLine(p1_out.GetX(fTelCS), p1_out.GetY(fTelCS),
                                         p2_out.GetX(fTelCS), p2_out.GetY(fTelCS));
          l_out->SetLineColor(1);
          l_out->SetLineStyle(2);

          l_in->Draw();
          l_out->Draw();
          objs->AddLast(l_in);
          objs->AddLast(l_out);

        }
      }
    }
  }

  // the camera support

  const Point pS(0. ,0. ,0. , fCameraSupportCoordinateSystem);

  const double dx = 5*cm;
  //double dy = 10*cm;
  const double length = 2*m;

  for (int i = 0; i < 2; ++i) {

    const double y1 = 30*cm;
    const double y2 = 30*cm + dx;
    const double z1 = 0*cm;
    const double z2 = 10*cm;
    const double x1 = 0*cm;
    const double x2 = length;

    TPolyLine3D* const up = new TPolyLine3D(4);
    TPolyLine3D* const down = new TPolyLine3D(4);
    up->SetLineColor(1);
    up->SetLineWidth(2);
    down->SetLineColor(1);
    down->SetLineWidth(2);

    double y = 0;
    double z = 0;
    for (int iy = 0; iy < 2; ++iy) {

      if (iy == 0)
        y = y1;
      else
        y = y2;

      if (i == 0)
        y *= -1;

      for (int iz = 0; iz < (do3D ? 2 : 1); ++iz) {

        if (iz == 0)
          z = z1;
        else
          z = z2;

        const Point p1(x1, y, z, fCameraSupportCoordinateSystem);
        const Point p2(x2, y, z, fCameraSupportCoordinateSystem);

        if (do3D) {  // DO THE 3D PLOT

          TPolyLine3D* const l = new TPolyLine3D(1);
          l->SetLineColor(1);
          l->SetLineWidth(2);
          l->SetPoint(0, p1.GetX(fTelCS), p1.GetY(fTelCS), p1.GetZ(fTelCS));
          l->SetPoint(1, p2.GetX(fTelCS), p2.GetY(fTelCS), p2.GetZ(fTelCS));
          l->Draw();
          objs->AddLast(l);

          up->SetPoint((1 - iy) * 2 + 1 - iz, p1.GetX(fTelCS), p1.GetY(fTelCS), p1.GetZ(fTelCS));
          down->SetPoint((1 - iy) * 2 + 1 - iz, p2.GetX(fTelCS), p2.GetY(fTelCS), p2.GetZ(fTelCS));

        } else {  // DO THE 2D PLOT

          TLine* const l = new TLine(p1.GetX(fTelCS), p1.GetY(fTelCS),
                                     p2.GetX(fTelCS), p2.GetY(fTelCS));
          l->SetLineColor(1);
          l->SetLineWidth(1);
          l->Draw();
          objs->AddLast(l);

          //up = new TLine();
          //down = new TLine();

        }

      }
    }

    if (do3D) {  // DO THE 3D PLOT
      up->Draw();
      objs->AddLast(up);
      down->Draw();
      objs->AddLast(down);
    } else {  // DO THE 2D PLOT
      //up->Draw();
      //down->Draw();
    }
  }

  return objs;
}