Modules/FdSimulation/TelescopeSimulatorKG/Camera.cc

Go to the documentation of this file.

//#define DEBUG_RT 1
#undef DEBUG_RT
#define DEBUG 1
//#undef 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 <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 <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 TelescopeSimulatorKG;
using namespace utl;


Camera::Camera(RandomEngine& rndm, const fdet::Telescope& tel) :
  fRandom(rndm),
  fTel(tel)
{
  const fdet::Camera& camera = tel.GetCamera();
  fTelCS = tel.GetTelescopeCoordinateSystem();
  fOrigin = tel.GetPosition();
  fRFocal = camera.GetRadiusFocal();
  fMercedesEff = camera.GetMercedesEfficiency();
  fSigma = camera.GetSigmaNormal();

  const double hmerc  = camera.GetMercedesHeight();
  const double ddmerc = camera.GetMercedesBase();   // half base length (full 9.2mm)

  // Warning: the values of SIN_?? may need some changes to simulate
  // correctly the camera (shadow and trace_surf).
  fSinAx = 0.2491262;
  fSinAy = 0.27966724;
  fSinBx = 0.25646528973;
  fSinBy = 0.26704475;

  // Camera
  SetMercedesParameter(hmerc, ddmerc);

  //------------ offset of the center of the camera ---------
  //fNx0 = 10.0;
  //fNy0 = 11.66666666666667;

  CalculateCameraSupportCoordinateSystem();

} // End of Camera::Camera


void
Camera::SetMercedesParameter(double height, double width)
{
  fR2 = 26.34*mm;

  // Find the position of the vertex in the reference frame of the
  // hexagon. See picture.

  // Full pixel vertex length
  //fR2 = sqrt (pixel209Pos.GetR (centerOfCameraCS));
  //cout << "SetMercParams : " << fR2/mm << endl;

  //----------------- Set mercedes parameters ------------------
  fHMerc = height;         //                    : height mercedes
  fD2 = fR2 * utl::kSqrt3; // 2 * cos(30) * r    : horizontal full pixel width
  fD1 = fD2 - 2.0*width;   //                    : horizontal pixel width (inside mercedes)
  fR1 = fD1/utl::kSqrt3;   //                    : hexagon vertex (inside mercedes)
  fNz = (fD2 - fD1)/2.0;   // == width           : mercedes width (!)
  fNt = fHMerc;            // == height
  const double aux = sqrt(fNz*fNz + fNt*fNt);
  fNz /= aux;              // (Nz, Nt) : normal on mercedes CHECK !! Seems to be wrong
  fNt /= aux;

  /*

        ^      /\
       n \    /  \
          \  /    \
           \/      \     height
           /        \
      ^   /          \
     z|  /            \
      |  --------------
      |       2xwidth
  <----
   t

 */
}// End of Camera::SetMercedesParameter


RTResult
Camera::Trace(const Photon& photonIn,
              Photon& photonOut,
              int& nreflections,
              int& col, int& row,
              double& cosTheta,
              bool doMercedes)
{
  nreflections = 0;

  //simulate camera, including/excluding the mercedes
  // cosTheta : photon hitting PMT angle
  // error codes:

  utl::Photon photonCamera;
  Vector normal;
  if (!TraceSurf(photonIn, photonCamera, normal)) {
#ifdef DEBUG
    if (TelescopeSimulatorKG::RayTracer::GetDebugLevel()>=1) 
      ERROR("RT::Camera> camera sphere not hit!");
#endif
    return eMissedFocalSurface;
  }
  
  Vector incidentPhotonDir(photonCamera.GetPosition() - fOrigin);
  incidentPhotonDir.Normalize();
  if (!FindPixelId(incidentPhotonDir, row, col)) {
    //the ray don't hit any pixel
    //cout << " no pixel hit " << endl;
    // but keep point on camera sphere anyway
    photonOut = photonCamera;
#ifdef DEBUG
    if (TelescopeSimulatorKG::RayTracer::GetDebugLevel()>=1) 
      INFO("RT::Camera> no pixle hit on camera!");
#endif
    return eMissedPixels;
  }

  // PixelId_new (pCamera-fOrigin, col, row); too slow ...

  // find the reference frame of PIXEL (oAux,mAux)
  //double aux1,aux2;
  const CoordinateSystemPtr& pxlCS = fTel.GetPixel(row, col).GetPixelCoordinateSystem();

  if (photonCamera.GetDirection().GetZ(pxlCS) >= 0) {
#ifdef DEBUG
    if (TelescopeSimulatorKG::RayTracer::GetDebugLevel()>=1) 
      INFO("RT::Camera> photon outgoing in pixel CS ???");
#endif
    return eBackscattered; // the light is going out from the pixel.
  }

  if (!doMercedes) {          //  no mercedes

    photonOut = photonCamera;
    cosTheta = -photonOut.GetDirection().GetCosTheta(pxlCS); // hitting the Camera

    return eOK; // hit camera after 0 reflections
  }

  // with mercedes

  // simulate mercedes
  utl::Photon photonPMT;
  const RTResult mercStatus = TraceMerc(pxlCS, photonCamera, photonPMT, nreflections);
  cosTheta = -photonPMT.GetDirection().GetCosTheta(pxlCS); // hitting the PMT
  
  if (mercStatus == eOK) { // HIT THE PMT after nreflections reflections
    
    double r_index_12 = 1.0/1.52;
    Vector normalPMT(0., 0., 1., pxlCS);
    RTFunctions::PhotonList photonsHittingPhotocathode;
    RTFunctions::Refraction(r_index_12, photonPMT, normalPMT, photonsHittingPhotocathode);    
    utl::Photon photonHittingPhotocathode = photonsHittingPhotocathode[0];
    
#ifdef DEBUG_RT
    int pixelId = (col-1)*fNRows + row;
    cout << "11111111 " << fTel.GetId() << " " << pixelId
         << " " << col
         << " " << row
         << " " << cosTheta
         << " " << photonHittingPhotocathode.GetWeight()
         << " " << photonPMT.GetWeight() << endl;
#endif

    photonOut = photonPMT;
    photonOut.SetWeight(photonHittingPhotocathode.GetWeight());

#ifdef DEBUG
    if (TelescopeSimulatorKG::RayTracer::GetDebugLevel()>=1) 
      {    ostringstream info;
    info << "RT::Camera>  HIT: photon: wIn=" << photonIn.GetWeight()
         << " wPMT=" << photonPMT.GetWeight()
         << " wSignal=" << photonOut.GetWeight()
         << " angle_cathode=" << acos(cosTheta)/deg << "deg "
         << " nrefl_merc=" << nreflections;
    INFO(info.str().c_str());
      }
#endif

    return eOK;
  }

#ifdef DEBUG
    if (TelescopeSimulatorKG::RayTracer::GetDebugLevel()>=1) 
      INFO("RT::Camera> nref<1 ! photon lost!");
#endif

  return mercStatus;  // lost on multi-reflection on mercedes (nreflections<0)
  
} // End of Camera::Trace

bool
Camera::FindPixelId(const Vector& incidentDir, int& row, int& col)
  const
{
  static const double fNx0 = 10.0;
  static const double fNy0 = 11.66666666666667; // optical axis
  static const double dOmega = fTel.GetCamera().GetEta();
  static const double dPhi = kSqrt3 * dOmega/2;
  static const int nCols = fTel.GetLastColumn() - fTel.GetFirstColumn() + 1; // 20
  static 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.0*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);
}


RTResult
Camera::TraceShadow(const utl::Photon& photonIn, const bool doSupport)
{
  // simulate shadow effect of the camera
  
  RTFunctions::IntersectionList intersections;
  if (RTFunctions::Sphere(fOrigin, fRFocal, photonIn, intersections)) {
    
    const Vector normal(intersections.back().second);

    const double x = normal.GetX(fTelCS);
    const double y = normal.GetY(fTelCS);

    if (-fSinAx < x && x < fSinBx&&
        -fSinAy < y && y < fSinBy) {

#ifdef DEBUG
    if (TelescopeSimulatorKG::RayTracer::GetDebugLevel()>=1) 
      INFO("RT::Camera> camera shadow!");
#endif

      return eShadowed; // the light hit the PMT camera
    }
  }
  
  Photon photonOut(intersections.back().first);

  if (doSupport) {

    Point pCamera(0., 0., 0., fCameraSupportCoordinateSystem);
    Vector nSupport(0., 0., 1., fCameraSupportCoordinateSystem);
    if (RTFunctions::Plane(pCamera, nSupport, photonIn, photonOut)) {

      Point pPlane = photonOut.GetPosition();

      double x = pPlane.GetX(fCameraSupportCoordinateSystem);
      double y = std::fabs(pPlane.GetY(fCameraSupportCoordinateSystem)); // symetric

      if (x < 0)
        return eOK;                       // camera support is only below the camera !

      if (y > 30.*cm && y < 35.*cm) {

#ifdef DEBUG
    if (TelescopeSimulatorKG::RayTracer::GetDebugLevel()>=1) 
        INFO("RT::Camera> camera support shadow!");
#endif

        return eShadowed;
      }
    }
  }

  return eOK;
} //End of Camera::TraceShadow


int
Camera::TraceSurf(const utl::Photon& photonIn, utl::Photon& photonOut,
                  Vector& /*normal*/)
{
  // find the point that the light hit the focal surface (pOut)
  // --------------
  //pIn = A point where the ray passes
  //nIn = direction of light propagation
  // ++++++++++++++y
  // pOut = Position that the ray hit the focal surface
  // nOut = ray direction.
  // **************
  // 0 -> The ray don't hit the camera.
  // 1 -> the ray hit the camera

  RTFunctions::IntersectionList intersections;
  if (RTFunctions::Sphere(fOrigin, fRFocal, photonIn, intersections) > 0) {
    
    photonOut = intersections.front().first;
    return 1;    // hit the camera sphere
    
    /*if(-fSinAx<nAux2[0] && nAux2[0]<fSinBx&&
      -fSinAy<nAux2[1] && nAux2[1]<fSinBy) {
    */
    //  }
  }
  
  return 0;  // miss the camera sphere
} // End of Camera::TraceSurf


RTResult
Camera::TraceMerc(const CoordinateSystemPtr& pxlCS,
                  const utl::Photon& photonIn, utl::Photon& photonOut,
                  int& nreflections)
{
  //simulate Mercedes
  //-------------
  // p = a point where the light passes
  // n = direction of incident light
  // ++++++++++++
  // ps = position at photo cathode surface.
  // ns = direction of light that hit photo cathode surface
  // ************
  // 0  -> the light don't hit the photo cathode
  // >0 -> number of times plus 1 that light hit mercedes bebfore hit
  //        photo cathode surface

  Vector dir(photonIn.GetDirection());

  Vector npl;
  Point o;
  utl::Photon point(photonIn); // starting point

  /*
  //---------- Verify if the light hit the considered pixel ---------
  ConstrPl (7, npl, o, pxlCS);                 // find the plane defined by the top of the mercedes.

  utl::Photon point;
  RTFunctions::Plane (o, npl, photonIn, point);  // move to plane above PMT (top of mercedes)

  //cout << " p " << point.GetX (fTelCS) << " " << point.GetY (fTelCS) << " " << point.GetZ (fTelCS) << endl;

  if (Hexagon(point.GetX(pxlCS), point.GetY(pxlCS),fR2) != eInside) {
    return(0); // the light is out of the limits of the pixel
  }
  */

  //------------- simulate refletions in the mercedes. -----------------
  nreflections = 0;
  int nloop = 0;
  int i_refl = -1;          // init
  utl::Photon point_next;        // the next point of ray-tracing
  int i = 0;
  do {                      // loop until the light hit photo cathode

    bool do_reflection = false;

    for (i = 0; i <= 6; ++i) { //find the plane that light hit.

      // 0 to photo cathode surface and 1,2,3,4,5 or 6 to mercedes

      if (!ConstructPlane(i, npl, o, pxlCS)) {
        ostringstream err;
        err << "This can never happen! "
            << "i: " << i
            << " nloop: " << nloop
            << " nrefl: " << nreflections
            << ". THROW EXCEPTION!";
        ERROR(err.str().c_str());
        throw(utl::AugerException(err.str()));
      }

      const double dist = RTFunctions::Plane(o, npl, point, point_next);
      const Point& pNext = point_next.GetPosition();

      /*
      cout << " dist: " << setw(8) << setprecision (3) << dist/mm
           << " i: " << i
           << " i_refl: " << i_refl
           << " nloop: " << nloop
           << " nrefl: " << nrefl
           << " side: " << Hexagon(pNext.GetX(pxlCS), pNext.GetY(pxlCS), fR1);
      */

      if (dist > 0 &&                                                   // photon will hit mercedes in the future
          Hexagon(pNext.GetX(pxlCS), pNext.GetY(pxlCS), fR1) == i &&    // intersection point is on right mercedes
          i != i_refl) {                                                // this is NOT the same mercedes we were just reflected by !
        do_reflection = true;
        //cout << " FOUND " << endl;
        break;
      }
      //cout << endl;

    } // for loop planes 0 - 6

    if (!do_reflection) {
#ifdef DEBUG
    if (TelescopeSimulatorKG::RayTracer::GetDebugLevel()>=1) 
      {      ostringstream err;
      err << " RT::Camera>  NO MERCEDES OR PMT HIT! PHOTON LOST! w=" << photonIn.GetWeight();
      INFO(err.str().c_str());
      }
#endif
      return eInvalidHitMercedes;
    }

    i_refl = i;                                    // remember the mercedes
    ++nloop;
    
    /*
      i==0: HIT THE PMT
    */

    if (i != 0) {                        // -> We hit a mercedes

      //RTFunctions::Normal(npl, fSigma, &fISeed, normal);
      const Vector normal = npl;
      const int status = RTFunctions::Reflection(point_next, normal, point_next);

      if (!status) {

#ifdef DEBUG
    if (TelescopeSimulatorKG::RayTracer::GetDebugLevel()>=1) 
      {        ostringstream err;
        err << "RT::Camera> get rid of this: photon hits mercedes from outside "
            << " " << npl*point_next.GetDirection() << " w=" << photonIn.GetWeight();
        INFO(err.str().c_str());
      }
#endif
        return eInvalidMercedes;
      }

      if (point_next.GetDirection().GetZ(pxlCS) >= 0) {         // backscattered from PMT?

#ifdef DEBUG
    if (TelescopeSimulatorKG::RayTracer::GetDebugLevel()>=1) 
      {        ostringstream err;
        err << "RT::Camera> backscattered from Mercedes? "
            << " theta="<< point_next.GetDirection().GetZ (pxlCS)
            << " w=" << photonIn.GetWeight();
        INFO (err.str().c_str());
      }
#endif
        return eBackscattered;                                  // this photon will not hit the PMT
      }

      point = point_next;
      point.SetWeight(point_next.GetWeight()*fMercedesEff);
      ++nreflections;                                                  // Do next reflection
    }

  } while (i != 0);                                             // WHILE we hit the PMT

  photonOut = point_next;
  
  return eOK;
  
} // End of Camera::TraceMerc


/*
int
Camera::PixelId_new(const Vector& vFocal, int& col, int& row)
{
  // find the pixel
  // --------------
  // vFocal: position at focal surface
  // ++++++++++++++
  // col,row -> identification of column and row of pixel
  // **************
  // 0 -> the position vFocal don't hit any pixel
  // 1 -> found a pixel

  //double omega, phi;
  double angleMin = 0;
  bool first = true;
  for (int iCol = 1; iCol <= fNCols; ++iCol) {
    for (int iRow = 1; iRow <= fNRows; ++iRow) {

      Vector pixelDir = - fTel.GetPixel(iRow, iCol).GetDirection();

      //double angle = Angle (pixelDir, vFocal);
      double angle = - std::fabs(pixelDir * vFocal);

      if (first) {

        angleMin = angle;
        col = iCol;
        row = iRow;
        first = false;

      } else {

        if (angle < angleMin) {
          angleMin = angle;
          col = iCol;
          row = iRow;
        }
      }
    }
  }

  const CoordinateSystemPtr& pxlCS = fTel.GetPixel(row, col).GetPixelCoordinateSystem();
  
  const Point pFocal = fOrigin + vFocal*fRFocal/vFocal.GetMag();
  
  if (eInside != Hexagon(pFocal.GetX(pxlCS), pFocal.GetY(pxlCS), fR2))
    return 0;

  return 1;
}
*/



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

   R: length of vertex


 */


Camera::HexagonDirection
Camera::Hexagon(double x, double y, 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

  double xl = std::fabs(x);
  double yl = std::fabs(y);
  double xm = utl::kSqrt3/2.0*r; // cos(30)*r
  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);
}


int
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 plano identified by id
  // 1 -> ok

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

  if (id > 0 && id < 7) {

    double dTheta = 60*deg;
    double theta = dTheta * id;
    double sinTheta = sin(theta); // y
    double cosTheta = cos(theta); // x

    n = Vector(-fNt*cosTheta, -fNt*sinTheta, fNz, pxlCS);
    p = Point(cosTheta*fD1/2., sinTheta*fD1/2., 0.0, pxlCS);

    return 1;
  }

  return 0;
}




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

  double xMin = 0, xMax = 0, yMin = 0, yMax = 0, zMin = 0, 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();

      double r_in = fR1;
      double r_out = fR2;

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

        double theta1 = 30*deg + 60*deg*i;
        double theta2 = 30*deg + 60*deg*(i+1);

        Point p1_in(cos(theta1) * r_in,
                    sin(theta1) * r_in,
                    0, pxlCS);
        Point p2_in(cos(theta2) * r_in,
                    sin(theta2) * r_in,
                    0, pxlCS);

        Point p1_out(cos(theta1) * r_out,
                     sin(theta1) * r_out,
                     fHMerc, pxlCS);
        Point p2_out(cos(theta2) * r_out,
                     sin(theta2) * r_out,
                     fHMerc, pxlCS);

        // pixel position in telCS
        double pxl_x = Point(0,0,0, pxlCS).GetX(fTelCS);
        double pxl_y = Point(0,0,0, pxlCS).GetY(fTelCS);
        double pxl_z = Point(0,0,0, pxlCS).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 *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(.01);
          p_pos->SetTextAngle(0);
          p_pos->Draw();
          objs->AddLast(p_pos);
        }

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

          TPolyLine3D *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 *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 *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 *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 *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);

        } // do 3D
      } // loop 6 mercedes
    } // loop row
  } // loop col

  // the camera support --------------------------------------------------------

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

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

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

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

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

    double y, z;
    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;

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

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

          TPolyLine3D *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 *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();

        }

      } // iz
    } // iy

    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();
    }
  } // iy

  return objs;
} // draw


void
Camera::CalculateCameraSupportCoordinateSystem()
{
  double theta = 16.*deg;
  double shift = fRFocal - 30.*cm;

  Vector vShift(0, 0, -1, fTelCS);
  CoordinateSystemPtr auxCS = CoordinateSystem::Translation(vShift*shift, fTelCS);
  fCameraSupportCoordinateSystem = CoordinateSystem::RotationY(theta, auxCS);
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "cd $AUGER_BASE && make"
// End: