/Lens.cc

Go to the documentation of this file.

#include "Lens.h"
#include "RTFunctions.h"

#include <TPolyLine3D.h>
#include <TObjArray.h>

#include <utl/Math.h>
#include <utl/MathConstants.h>
#include <utl/AugerUnits.h>

#include <utl/Photon.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/CoordinateSystemPtr.h>

#include <utl/RandomEngine.h>
#include <CLHEP/Random/Randomize.h>

#include <fdet/Telescope.h>
#include <fdet/Corrector.h>

#include <det/Detector.h>

#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>

#include <iostream>
#include <sstream>
#include <utility>

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

using CLHEP::RandFlat;


Lens::Lens(utl::RandomEngine& rndm,
           const fdet::Telescope& detTel,
           const bool simulateHaloEffects,
           const double lensIncreaseRefl,
           const double lensPosition,
           const double minLensThickness,
           const double torusRadius,
           const double tubeRadius,
           const double torusZ0) :
  fRandom(rndm),
  fSimulateHaloEffects(simulateHaloEffects),
  fIncreaseReflection(lensIncreaseRefl),
  fPosZ(lensPosition),
  fMinLensThickness(minLensThickness),
  fTorusRadius(torusRadius),
  fTubeRadius(tubeRadius),
  fTorusZ0(torusZ0),
  fTelCS(detTel.GetTelescopeCoordinateSystem()),
  fOrigin(0, 0, fPosZ, fTelCS),
  fTorusOrigin(0, 0, fTorusZ0, fTelCS),
  fPosPlane(0, 0, fPosZ, fTelCS),
  fR1(detTel.GetCorrector().GetInnerRadius()),
  fR2(detTel.GetCorrector().GetOuterRadius()),
  fRefractiveIndex(detTel.GetCorrector().GetRefractiveIndex()),
  fTransmittance(detTel.GetCorrector().GetTransmittance()),
  fSigRho(detTel.GetCorrector().GetSigmaNormal()),
  fNbug(0)
{
  // lens defines the entrance of the telescope
  // center of curvature of camera and mirror are here
  // 0 is at the edge of the lens elements at 850mm (the one closer to the filter)

  /*fPosZ = -0.*cm;
  fMinLensThickness = 6.628*mm;  // Julias first calculation
  fTorusRadius = 795.*mm;  // FD paper
  fTubeRadius = 8380.*mm;  // FD paper
  fTorusZ0 = fPosZ + ftubeRadius - 0.18*mm;  // Julias first calculation

  fLXTorusRadius = 795.27*mm;  // value used in LX module (spherical)
  fLXTubeRadius = 8383.*mm;  // value used in LX module (spherical)
  fLXz0 = 8380.9*mm;

  fMinLensThickness = 6.13822*mm;//6.5*mm;
  fTorusRadius = 795.27*mm;//795.*mm;//795.*mm;//849.*mm;//795.*mm;//795.*mm;//850.*mm;//804.10*mm;
  fTubeRadius = 8383*mm;//9002.67*mm;//8832.98*mm;//10169.7*mm;//21436.5*mm;//14524.6*mm;//6130.*mm;//ftubeRadius;
  const double tmp_torusZ0 = ftubeRadius + 5.95956*mm;//9002.5*mm;//8826.25*mm;//10162.8*mm;//21429.7*mm;//14517.8*mm ;//ftubeRadius + 6.9*mm;//8386.775*mm;
  ftorusZ0 = fPosZ + tmp_torusZ0;*/

  //fTorusOrigin = Point(0, 0, fPosZ + ftubeRadius - 0.18*mm, fTelCS);
  //fPosPlane = Point(0, 0, fPosZ - fminLensThickness, fTelCS);
  //fPosPlane = Point(0, 0, -fminLensThickness, fTelCS);
  // const double Tm = corrector.GetMeanLensThickness();      // average lens thickness

  // air refractive index
/*#warning this is inconsistent throughout the raytracing! CHECK.
#warning Here refractive index of air is taken from the time-dependent Atmosphere
  const atm::ProfileResult& refractiveIndexVsHeight =
    det::Detector::GetInstance().GetAtmosphere().EvaluateRefractionIndexVsHeight();
  fIndexAir = refractiveIndexVsHeight.Y(refractiveIndexVsHeight.MinX());*/
}


#if 0
RTNext
Lens::Trace(const utl::Photon& photonIn, utl::Photon& photonOut)
{
  cout << "tracing without torus!" << endl;

  // some helper variables
  const Vector nPlane(0, 0, 1, fTelCS);
  const double lambda = photonIn.GetWavelength();
  const double r_index = fRefractiveIndex.Y(lambda) / fIndexAir;
  double weight = photonIn.GetWeight();  // input weight, this will get reduced via absorption

  // propagate photon to lens "plane" -> photonLens
  // at this stage the lens is just a plane at fPosZ distance from aperture
  // fOrigin is Point(0, 0, fPosZ, fTelCS)
  utl::Photon photonLens;
  RTFunctions::Plane(fPosPlane, nPlane, photonIn, photonLens);

  bool inLens = false;  // photon not yet in lens material

  /* the corrector lens gaps are implemented, but there are no vertical optical boundaries implemented.
    Only the front and back side of the lens plus the inner/outer radius + gaps.
    Effectively: photons that travel beyond R2 are absorbed, photons that travel beyond R1 or hit one of the gaps,
    pass out of the lens back into the telescope. */

  //double travelPath = 0;  // sum up the distance traveled by the photon in the glass

  int nDist = 0;
  int count = 0;

  do {  // internal reflections in the glass

    ++count;
    const Vector& nDir = photonLens.GetDirection();
    const bool backwards = nDir.GetZ(fTelCS) > 0;  // towards the aperture !!

    if (!fSimulateHaloEffects && backwards) {
      // the light is going out of the telescope
      photonOut = photonLens;
      return eLost;
    }

    const Point& pLens = photonLens.GetPosition();
    const double x = pLens.GetX(fTelCS);
    const double y = pLens.GetY(fTelCS);
    const double r2 = x*x + y*y;
    const double r = sqrt(r2);

    /*if (r > fR1 && r < fR2) {
      photonOut = photonLens;
      return eLost;  // mask the lens
    }*/
    if (r > fR2) {
      photonOut = photonLens;
      return eLost;  // the light doesn't pass through the aperture window
    }

    if (r < fR1) {
      photonOut = photonLens;
      if (inLens) {
        //cout << "first inLens, now in aperture - bug! for now -> return eLost" << endl;
        ++fNbug;
        return eLost;
      }
      // the light went through the aperture window,
      return backwards ? eFilter : eCameraShadow;
      // to mask the center of the aperture use instead:
      //return eLost;
    }

    // corrector ring gaps
    const double phiAngle = atan(abs(y/x));
    if (fSimulateHaloEffects && phiAngle - floor(phiAngle / (kPi/12)) * (kPi/12) < 0.265*deg /*gap*/) {
      photonOut = photonLens;
      // avoid photons going from lens glass directly into the gap (would result in wrong direction)
      if (inLens) {
        //cout << "first inLens, now in gap - bug! for now -> return eLost" << endl;
        ++fNbug;
        return eLost;
      }
      return backwards ? eFilter : eCameraShadow;
    }

    const Vector normal = backwards ?
      (inLens ? CurvaturePrototype(x,y) : nPlane) :
      (inLens ? nPlane : CurvaturePrototype(x, y));
    cout << " we are doing this: curvaturePrototype can be used!" << endl;

    // calculate refraction at optical boundary
    RTFunctions::PhotonList photonsRefracted;
    RTFunctions::Refraction(inLens ? r_index : 1/r_index, photonLens, normal, photonsRefracted);

    if (photonsRefracted.size() == 1) {  // total internal reflection

      if (!fSimulateHaloEffects) {
        photonOut = photonsRefracted[0];
        return eLost;
      }

      if (!inLens) {
        cout << "total internal reflection but !inLens ---------- this can not happen!" << endl;
        photonOut = photonsRefracted[0];
        return backwards ? eCameraShadow : eFilter;
      }

      // photon is total internal reflected in lens

      photonLens = photonsRefracted[0];

      // -> propagate to next optical surface

    } else {  // partial transmission/reflection

      // modify reflectivity (ad-hoc) ONLY when coming from the outside! Internal reflections are not affected.
      if (!inLens && fIncreaseReflection != 1) {
        const double weightRefl = photonsRefracted[1].GetWeight();
        const double inWeight = weight;
        const double newRefl = min(inWeight, weightRefl * fIncreaseReflection);
        const double newRefr = inWeight - newRefl;
        photonsRefracted[1].SetWeight(newRefl);
        photonsRefracted[0].SetWeight(newRefr);
      }
      // end modify reflectivity --

      if (fSimulateHaloEffects) {  // in the halo case we do MC raytracing of one photon

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

        if (tmpRandNr > photonsRefracted[0].GetWeight() / weight) {

          // -> reflection
          if (!inLens) {  // photon does not enter the lens
            photonOut = photonsRefracted[1];
            photonOut.SetWeight(weight);
            //cout << "photon does not enter the lens, count = " << count << endl;
            return backwards ? eCameraShadow : eFilter;
          }

          photonLens = photonsRefracted[1];
          photonLens.SetWeight(weight);

          // -> we are internal in lens glass

          // -> propagate to next optical interface ...

        } else {  // end reflection, start transmission

          if (!inLens) {

            inLens = true;  // now we are in the lens !

            photonLens = photonsRefracted[0];
            photonLens.SetWeight(weight);
            //cout << "now in lens, count: " << count << endl;

            // -> propagate to next optical interface

          } else {  // now we are leaving the lens !

            photonOut = photonsRefracted[0];
            photonOut.SetWeight(weight);
            //cout << "now leaving lens, count: " << count << endl;
            // debug
            //cout << "x: " << x << ", y: " << y << endl;
            return backwards ? eFilter : eCameraShadow;

          }

        }

      } else {  // end do Halo, start non-halo

        // in the non-halo mode we are only considering the weight-loss of the photon moving in forward direction!

        // only (forward) transmission

        if (inLens) {  // photon is leaving the lens
          photonOut = photonsRefracted[0];
          return eCameraShadow;
        }

        inLens = true;
        photonLens = photonsRefracted[0];
        weight = photonLens.GetWeight();  // update weight

        // -> propagate to next optical interface ...

      }

    }

    // propagate photon to next optical interface

    const Vector& nDirLens = photonLens.GetDirection();
    const bool backwardsLens = nDirLens.GetZ(fTelCS) > 0;  // towards the aperture !!
    const double cosThetaPhoton = (backwardsLens ? +1 : -1) * nDirLens.GetCosTheta(fTelCS);
    const double z = ProfilePrototype(r);  // glass thickness
    const double distance = z / cosThetaPhoton;  // thickness of the glass along nDirLens

    // RU: save track point for drawing
    if (fPhotonTrack)
      fPhotonTrack->AddPoint(pLens);

    // calculate the position of light at the second surface of the lens.
    Point pNext = pLens + distance * nDirLens;
    if (!backwardsLens)  // in this case we are now at the plane inner surface of the lens at fPosZ from the filter
      pNext = Point(pNext.GetX(fTelCS), pNext.GetY(fTelCS), fPosZ, fTelCS);

    photonLens.SetPosition(pNext);

    weight *= exp(distance / (1*cm) * log(fTransmittance.Y(lambda)));
    photonLens.SetWeight(weight);

    ++nDist;
    //cout << " -> distance=" << distance/mm << " mm nDist=" << nDist << endl;

    // mark photons and plot them later
    if (nDist == 3 /*nDist >= 23*/) {
      const int everyOnehundred = RandFlat::shoot(&fRandom.GetEngine(), 0, 100);
      if (everyOnehundred == 4)
        photonLens.SetSource(3333);  // secret number
    }

  } while (weight > 1e-20);

  // too many multiple reflections -> absorbed
  photonOut = photonLens;
  return eLost;
}
#endif 

RTNext
Lens::TraceWithTorus(const utl::Photon& photonIn, utl::Photon& photonOut)
{
  const bool printLens = false;
  // some helper variables
  const double lambda = photonIn.GetWavelength();
  const double r_index = fRefractiveIndex.Y(lambda) / fIndexAir;
  //const double r_index = 1.2;
  //cout << "refractive index r_index = " << r_index << endl;
  double weight = photonIn.GetWeight();  // input weight, this will get reduced via absorption

  const Vector inwards(0, 0, -1, fTelCS);
  const bool outwards  = photonIn.GetDirection().GetZ(fTelCS) > 0;
  const Vector nPlane(0, 0, 1, fTelCS);

  // propagate photon to lens surface -> photonLens
  // fOrigin is Point(0, 0, fPosZ, fTelCS)

  utl::Photon photonLens;
  if (outwards)
    RTFunctions::Plane(fPosPlane, nPlane, photonIn, photonLens);
  else {
    RTFunctions::IntersectionList intersections;
    RTFunctions::Torus(fTorusOrigin, inwards, fTorusRadius, fTubeRadius, photonIn, intersections);
    if (intersections.empty()) {
      ERROR("no intersection with torus???");
      return eLost;
    }

    unsigned int indexOfSmallestZ = 4;  // maximum number of intersections is 4 (0 to 3)
    // curved lens surface has the smallest z values on the torus
    ChoosePhysicalIntersection(indexOfSmallestZ, intersections);

    if (indexOfSmallestZ == 4) {
      ERROR("smallestZ was not found -> indexOfSmallestZ is not correct!");
      return eLost;
    }

    photonLens = photonIn;
    photonLens.SetPosition(intersections[indexOfSmallestZ].first.GetPosition());
  }

  bool inLens = false;  // photon not yet in lens material

  /* the corrector lens gaps are implemented, but there are no vertical optical boundaries implemented.
    Only the front and back side of the lens plus the inner/outer radius + gaps.
    Effectively: photons that travel beyond R2 are absorbed, photons that travel beyond R1 or hit one of the gaps,
    pass out of the lens back into the telescope. */

  //double travelPath = 0;  // sum up the distance traveled by the photon in the glass

  int nDist = 0;
  int count = 0;

  do {  // internal reflections in the glass

    ++count;
//    const Vector& nDir = photonLens.GetDirection();
//    const bool backwards = nDir.GetZ(fTelCS) > 0;  // towards the aperture !!
    const bool outwards = photonLens.GetDirection().GetZ(fTelCS) > 0;// towards the aperture !!

    if (!fSimulateHaloEffects && outwards) {
      // the light is going out of the telescope
      photonOut = photonLens;
      return eLost;
    }

    const Point& pLens = photonLens.GetPosition();
    const double pos_x = pLens.GetX(fTelCS);
    const double pos_y = pLens.GetY(fTelCS);
    //const double pos_z = pLens.GetZ(fTelCS);
    const double pos_r2 = Sqr(pos_x) + Sqr(pos_y);
    const double pos_r = sqrt(pos_r2);

    /*if (pos_r > fR1 && pos_r < fR2) {
      photonOut = photonLens;
      return eLost;     //mask the lens
    }*/

    if (pos_r > fR2) {
      photonOut = photonLens;
      return eLost;  // the light doesn't pass through the aperture window
    }

    if (pos_r < fR1) {
      photonOut = photonLens;
      if (inLens) {
        //cout << "first inLens, now in aperture - bug! for now -> return eLost" << endl;
        ++fNbug;
        return eLost;
      }
      if (!outwards && printLens)
        cout << 0 << ',';
      // the light went through the aperture window,
      return outwards ? eFilter : eCameraShadow;

      // to mask the center of the aperture use instead:
      // return eLost;
    }

    const bool debug = false;

    // corrector ring gaps, see FD NIM paper for numbers (Fig. 11)
    const double ringElementAngle = kPi / 12; // = 15 deg
    const double phiAngle = atan2(abs(pos_y), abs(pos_x)) + (ringElementAngle / 2);
    const double gapWidth = ((2 * kPi * 850*mm) - (24 * 218*mm)) / 24; // = 4.5mm, fig.11 FD NIM
    //const double gapWidth = 19 * mm; // width of the tape on the lens gaps, see FD Task Log 2011-11-28
    const double gapAngularWidth = sin(gapWidth / (975*mm)); // 975mm = middle of corrector ring
    //if (fSimulateHaloEffects && abs(phiAngle - round(phiAngle / ringElementAngle) * ringElementAngle) < (gapAngularWidth / 2)) {
    if (abs(phiAngle - round(phiAngle / ringElementAngle) * ringElementAngle) < (gapAngularWidth / 2)) {
      if (debug) {
        cout << pLens.GetX(fTelCS) / mm << ','
             << pLens.GetY(fTelCS) / mm << ','
             << pLens.GetZ(fTelCS) / mm << '\n';
      }
      photonOut = photonLens;
      //return eLost;  // make gaps in the Lens absorbing 
      // avoid photons going from lens glass directly into the gap (would result in wrong direction)
      if (inLens) {
        //cout << "first inLens, now in gap - bug! for now -> return eLost" << endl;
        ++fNbug;
        return eLost;
      }
      //ERROR("hit the corrector gap");
      if (!outwards && printLens)
        cout << 2 << ',';
      return outwards ? eFilter : eCameraShadow;
    }

    const Vector nTorus = RTFunctions::TorusNormal(fTorusOrigin, inwards, pLens, fTorusRadius, fTubeRadius);
    // for test purposes use "old" function instead
    //const Vector nTorus = CurvaturePrototype(pos_x, pos_y);

    // depending on the incoming direction take normal on the plane or curved side of the lens
    const Vector normal = outwards ? (inLens ? nTorus : nPlane) : (inLens ? nPlane : nTorus);

    // debug
    /*cout << "backwards = " << backwards << ", inLens = " << inLens << ", normal = (" << normal.GetX(fTelCS) << ","
         << normal.GetY(fTelCS) << ","
         << normal.GetZ(fTelCS) << ") in point (" << pos_x/mm << "," << pos_y/mm << "," << pos_z/mm << ") "
            "at r = " << pos_r/mm << ", old z value = " << ProfilePrototype(pos_r)/mm << endl;
    cout << "for root " << backwards << "\t" << inLens << "\t" << normal.GetX(fTelCS) << "\t"
         << normal.GetY(fTelCS) << "\t"
         << normal.GetZ(fTelCS) << "\t" << pos_x/mm << "\t" << pos_y/mm << "\t" << pos_z/mm << "\t"
         << pos_r/mm << "\t" << ProfilePrototype(pos_r)/mm << endl;*/

    // calculate refraction at optical boundary
    RTFunctions::PhotonList photonsRefracted;
    RTFunctions::Refraction(inLens ? r_index : 1 / r_index, photonLens, normal, photonsRefracted);

    if (photonsRefracted.size() == 1) {  // total internal reflection

      if (!fSimulateHaloEffects) {
        photonOut = photonsRefracted[0];
        return eLost;
      }

      if (!inLens) {
        ERROR("total internal reflection but !inLens ---------- this can not happen!");
        photonOut = photonsRefracted[0];
        //return outwards ? eCameraShadow : eFilter;
        return eLost;
      }

      // photon is total internal reflected in lens

      photonLens = photonsRefracted[0];

      // -> propagate to next optical surface

    } else {  // partial transmission/reflection

      // modify reflectivity (ad-hoc) ONLY when coming from the outside! Internal reflections are not affected.
      if (!inLens && fIncreaseReflection != 1) {
        const double weightRefl = photonsRefracted[1].GetWeight();
        const double inWeight = weight;
        const double newRefl = min(inWeight, weightRefl * fIncreaseReflection);
        const double newRefr = inWeight - newRefl;
        photonsRefracted[1].SetWeight(newRefl);
        photonsRefracted[0].SetWeight(newRefr);
      }
      // end modify reflectivity --

      if (fSimulateHaloEffects) {  // in the halo case we do MC raytracing of one photon
        const double tmpRandNr = RandFlat::shoot(&fRandom.GetEngine(), 0, 1);
        if (tmpRandNr > photonsRefracted[0].GetWeight() / weight) {
          // -> reflection
          if (!inLens) {  // photon does not enter the lens
            photonOut = photonsRefracted[1];
            photonOut.SetWeight(weight);
            //cout << "photon does not enter the lens, count = " << count << endl;
            /*// mask lens
            if (backwards)
              return eLost;
            else*/

            // diffusive component
            const double amountOfStrayLight = 0;  // 0% of the light will go in "lambertian directions"
            if (amountOfStrayLight > 0) {
              double maxTheta = 25*deg;
              const Vector lambertdiffusion =
                RTFunctions::LambertDiffusion(fRandom , photonOut.GetDirection(), normal, amountOfStrayLight, maxTheta);
              photonOut.SetDirection(lambertdiffusion);
            }
            if (!outwards && printLens)
              cout << 1 << ',';
            return outwards ? eCameraShadow : eFilter;
          }

          photonLens = photonsRefracted[1];
          photonLens.SetWeight(weight);

          // -> we are internal in lens glass
          // -> propagate to next optical interface ...

        } else {  // end reflection, start transmission

          if (!inLens) {

            inLens = true;  // now we are in the lens !

            photonLens = photonsRefracted[0];
            photonLens.SetWeight(weight);
            //cout << "now in lens, count: " << count << endl;

            // -> propagate to next optical interface

          } else {  // now we are leaving the lens !

            photonOut = photonsRefracted[0];
            photonOut.SetWeight(weight);
            //cout << "now leaving lens, count: " << count << endl;
            /*// mask lens
            if (!backwards)
              return eLost;
            else*/

            // diffusive component
            //const double amountOfStrayLight = 0.05;  // 5% of the light will go in "lambertian directions"
            const double amountOfStrayLight = 0;  // 0% of the light will go in "lambertian directions"
            double maxTheta = 25*deg;
            const Vector lambertdiffusion =
              RTFunctions::LambertDiffusion(fRandom, photonOut.GetDirection(), normal, amountOfStrayLight, maxTheta);
            photonOut.SetDirection(lambertdiffusion);
 
            if (!outwards && printLens)
              cout << 1 << ',';
            return outwards ? eFilter : eCameraShadow;
          }
        }
      } else {  // end do Halo, start non-halo

        // in the non-halo mode we are only considering the weight-loss of the photon moving in forward direction!
        // only (forward) transmission

        if (inLens) {  // photon is leaving the lens
          photonOut = photonsRefracted[0];
          return eCameraShadow;
        }

        inLens = true;
        photonLens = photonsRefracted[0];
        weight = photonLens.GetWeight();  // update weight

        // -> propagate to next optical interface ...
      }
    }

    Point pNext = pLens;
    const bool outwardsLens = photonLens.GetDirection().GetZ(fTelCS) > 0;  // towards the aperture !!
    if (inLens) {
      if (outwardsLens) {  // intersection with torus
        RTFunctions::IntersectionList intersections;
        RTFunctions::Torus(fTorusOrigin, inwards, fTorusRadius, fTubeRadius, photonLens, intersections);

        if (intersections.empty()) {
          ERROR("inside the lens and flying outwards - but no intersection with torus???");
          return eLost;
        }

        unsigned int indexOfSmallestZ = 4;  // maximum number of intersections is 4 (0 to 3)
        ChoosePhysicalIntersection(indexOfSmallestZ, intersections);

        if (indexOfSmallestZ == 4) {
          ERROR("smallestZ was not found -> indexOfSmallestZ is not correct!");
          return eLost;
        }
        // curved lens surface has the smallest z values on the torus
        pNext = intersections[indexOfSmallestZ].first.GetPosition();

      } else { // intersection with plane
        Photon tmpPhoton = photonLens;
        RTFunctions::Plane(fPosPlane, nPlane, photonLens, tmpPhoton);
        pNext = tmpPhoton.GetPosition();

      }
      photonLens.SetPosition(pNext);
    } else {
      ERROR("propagate to next optical interface (inside lens) but photon is not inLens!!!!!");
      return eLost;
    }

    const Vector& pathInLens = pNext - pLens;
    const double distance = pathInLens.GetMag();
    weight *= exp(distance / (1*cm) * log(fTransmittance.Y(lambda)));

    ++nDist;
    if (nDist > 100) {
      ERROR("so many reflections? CHECK?!?!");
      return eLost;
    }

    // mark photons and plot them later
    if (nDist == 3 /*nDist >= 23*/) {
      const int everyOnehundred = RandFlat::shoot(&fRandom.GetEngine(), 0, 100);
      if (everyOnehundred == 4)
        photonLens.SetSource(3333);  // secret number
    }

  } while (weight > 1e-20);

  ERROR("Lost photon due to too many multiple reflections in the lens");
  // too many multiple reflections -> absorbed

  photonOut = photonLens;
  return eLost;
}


TObjArray*
Lens::Draw()
{
  TObjArray* const objs = new TObjArray();

  const int color = 1;
  const int width = 2;

  const int nSeg = 15;
  for (int i = 0; i < nSeg; ++i) {

    const double theta_1 = 360*deg / nSeg * i;
    const double theta_2 = 360*deg / nSeg * (i + 1);

    TPolyLine3D* const l_in = new TPolyLine3D(2);
    l_in->SetPoint(0, cos(theta_1)*fR1, sin(theta_1)*fR1, fPosZ);
    l_in->SetPoint(1, cos(theta_2)*fR1, sin(theta_2)*fR1, fPosZ);
    l_in->SetLineColor(color);
    l_in->SetLineWidth(width);
    l_in->Draw();
    objs->AddLast(l_in);

    TPolyLine3D* const l_out = new TPolyLine3D(2);
    l_out->SetPoint(0, cos(theta_1)*fR2, sin(theta_1)*fR2, fPosZ);
    l_out->SetPoint(1, cos(theta_2)*fR2, sin(theta_2)*fR2, fPosZ);
    l_out->SetLineColor(color);
    l_out->SetLineWidth(width);
    l_out->Draw();
    objs->AddLast(l_out);

    TPolyLine3D* const l_seg = new TPolyLine3D(2);
    l_seg->SetPoint(0, cos(theta_1)*fR1, sin(theta_1)*fR1, fPosZ);
    l_seg->SetPoint(1, cos(theta_1)*fR2, sin(theta_1)*fR2, fPosZ);
    l_seg->SetLineColor(color);
    l_seg->SetLineWidth(width);
    l_seg->Draw();
    objs->AddLast(l_seg);
  }
  return objs;
}


/*
 new KG default from GAP 1999-12,1999-13, 1999-14 (older prototypes)
 GAP 2000-003 describes also the final lens design !!!!!!!!!!!!!!!!!

 the function CurvaturePrototype and ProfilePrototype
 are derived from Tilo Waldenmaier, not sure how this relates to reality...

 The real lenses are torus shapes as explained in gap 2000-003 and in the FD-Paper NIM 2012
*/


// thikness profile as a function of radial distance
Vector
Lens::CurvaturePrototype(const double x, const double y)
  const
{
  const double r = sqrt(x*x + y*y);

  static const double f = 3.4*m - 1.743*m;
  static const double Rd = 0.85*m;
  static const double A = 3/2. * Rd * Rd;
  static const double n = 1.5;

  static const double denominator = 32 * (n - 1) * f * f * f;
  static const double a1 = 1 / denominator;  // r^4
  static const double a2 = A / denominator;  // r^2

  const double dTdr = 4 * a1 * r * r * r - 2 * a2 * r;
  //dTdr *= 0.9;  // this would be a bit better

  // Tilos KA prototyp
  //const double dTdr = 6 * 4e-18*mm * pow(r, 5) * pow(1*mm, -6);

  // new KG default from GAP ...

  double nr = dTdr;
  double nz = 1;
  const double norm = sqrt(nr*nr + nz*nz);

  nr /= norm;
  nz /= norm;

  const double phi = atan2(y, x);

  return Vector(-nr * cos(phi), -nr * sin(phi), nz, fTelCS);

  // old FDSim code:
  /*const double r2 = r*r;
  const double aux = 2*(-1.2461221e-02 + 2 * 8.6476346e-03 * r2 + 3*2.0361284e-03 * r2*r2);

  const Vector norm = Normalized(Vector(-x*aux, -y*aux, 1.0, fTelCS));

  return norm;*/
}


void
Lens::ChoosePhysicalIntersection(unsigned int& index,
                                 const RTFunctions::IntersectionList& intersections)
  const
{
  double smallestZ = 10 * fTubeRadius;  // definitely outside the torus and outside the telescope
  for (unsigned int i = 0, n = intersections.size(); i < n; ++i) {
    const double tmp_z = intersections[i].first.GetPosition().GetZ(fTelCS);
    if (tmp_z < smallestZ) {
      smallestZ = tmp_z;
      index = i;
    }
  }
}

#if 0
// thikness profile as a function of radial distance
double
Lens::ProfilePrototype(const double r)
  const
{
  const double r2 = r * r;

  // Tilos KA prototyp
  /*r2 /= (mm*mm);
  const double z = 4e-18*mm * pow(r2, 3) + 0.91*mm;*/

  static const double f = 3.4*m - 1.743*m;
  static const double Rd = 0.85*m;
  static const double A = 3/2. * Rd * Rd;
  static const double n = 1.5;

  static const double denominator = 32 * (n - 1) * f * f * f;
  static const double a1 = 1 / denominator;  // r^4
  static const double a2 = A / denominator;  // r^2

  static const double z0 = 10*mm;

  const double z = z0 + a1 * r2 * r2 - a2 * r2;

  // old FDSim code:
  /*const double aux = 2*(-1.2461221e-02 + 2*8.6476346e-03*r2 + 3*2.0361284e-03*r2*r2);

  const Vector norm = Normalized(Vector(-x*aux, -y*aux, 1.0, fTelCS));

  const double z = sqrt(pow(norm.GetX(fTelCS), 2) + pow(norm.GetY(fTelCS), 2));*/

  return z;
}
#endif


Lens::~Lens()
{
#if 0
  if (fNbug) {
    //INFO("#######################################################");
    ostringstream bugMsg;
    bugMsg << fNbug << " photons lost because of \'first inLens, now in aperture\'-bug";
    INFO(bugMsg.str());
    //INFO("#######################################################");
  }
#endif
}