Modules/FdSimulation/TelescopeSimulatorKG2/Mirror.cc

Go to the documentation of this file.

#include "Mirror.h"
#include "RTFunctions.h"

#include <utl/Math.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/RandomEngine.h>
#include <utl/TabulatedFunction.h>

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

#include <CLHEP/Random/Randomize.h>

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

#include <iostream>

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

using CLHEP::RandFlat;
using CLHEP::RandGauss;


// these are constants used for the segmentation of the mirror sphere
const int gNSegments = 6;
const double gMirrorWidth = 3.5*m;


Mirror::Mirror(utl::RandomEngine& rndm,
               const fdet::Telescope& tel,
               const double mirrorSize,
               const double mirrorSegmentSigma,
               const double mirrorRadiusSigma,
               const double mirrorAbsorptionTop,
               const double mirrorAbsorptionBot,
               const utl::TabulatedFunction* const mirrorDiffusionTop,
               const utl::TabulatedFunction* const mirrorDiffusionBot) :
  fRandom(rndm),
  fMirrorSize(mirrorSize),
  fMirrorSegmentSigma(mirrorSegmentSigma),
  fMirrorRadiusSigma(mirrorRadiusSigma),
  fMirrorAbsorptionTop(mirrorAbsorptionTop),
  fMirrorAbsorptionBot(mirrorAbsorptionBot),
  fMirrorDiffusionTop(mirrorDiffusionTop),
  fMirrorDiffusionBot(mirrorDiffusionBot),
  fTelCS(tel.GetTelescopeCoordinateSystem()),
  fOrigin(0,0,0, fTelCS),
  fReflectivity(tel.GetMirror().GetReflectivity()),
  fRCurv(tel.GetMirror().GetRadiusOfCurvature()),
  fSigma(tel.GetMirror().GetSigmaNormal())
{ }


RTNext
Mirror::Trace(const utl::Photon& photonIn, utl::Photon& photonOut)
{
  if (photonIn.GetDirection().GetZ(fTelCS) > 0) {
    photonOut = photonIn;
    ostringstream err;
    err << "Missed Mirror : nDir>0";
    ERROR(err);
    return eLost;
  }

  RTFunctions::IntersectionList intersections;
  const int sphereStatus = RTFunctions::Sphere(fOrigin, fRCurv, photonIn, intersections);
  if (sphereStatus <= 0) {
    photonOut = photonIn;
    ostringstream err;
    err << "Missed Mirror : " << sphereStatus << " !!!";
    ERROR(err);
    return eLost;
  }

  // always choose latter intersection, we want to hit inside of sphere
  // photonOut is the photon when it hits the mirror
  photonOut = intersections.back().first;
  Vector normal = intersections.back().second;

  // dimensions of the mirror +/- 33.07 deg in latitude and longitude by default (change in TelescopeSimulatorKG2 --> MirrorSize)
  // or use GetMirrorSegment(posOnSphere) that was implemented for the deformed mirror
  const Point posOnSphere = photonOut.GetPosition();

  //position of photon on mirror
  //latitude --> used to describe changing dust layer across mirror in vertical direction
  const double horizontalPosOnMirror = asin(posOnSphere.GetY(fTelCS) / fRCurv);  // longitude
  const double verticalPosOnMirror = asin(posOnSphere.GetX(fTelCS) / fRCurv);  // latitude

  // if mirror size is set to -1 in TelescopeSimulatorKG2.xml.in --> no mirror boundary, i.e. mirror is a complete hemisphere
  // if a physical mirror size is given (0 deg < mirrorSize < 90 deg), check now if photon hits mirror
  //if (fMirrorSize > 0) {
  if (false) {
    if (abs(horizontalPosOnMirror) > fMirrorSize || abs(verticalPosOnMirror) > fMirrorSize) {
      //cout << posOnSphere.GetX(fTelCS) / mm << ','
      //     << posOnSphere.GetY(fTelCS) / mm << ','
      //     << posOnSphere.GetZ(fTelCS) / mm << '\n';
      return eLost;
    }
  }
  
  if (true) { // mirror gaps
    bool lost = false;
    const bool HARDDEBUG = false;
    if (1) {
      // gaps between mirror segments for KA mirrors
      const int nSeg = 3;  // number of segments per row/col per Quadrant 
      const double mirrorGapWidth =  5*mm;  // educated guess
      const double angGapWidth = 2 * fRCurv * asin(0.5*mirrorGapWidth / fRCurv);  // approx 10 deg at equator
      const double mirrorSegmentHeight = 650*mm;  // = width of segment at equator
      const double mirrorSegmentAngWidth = 2 * fRCurv * asin(0.5*mirrorSegmentHeight / fRCurv);

      /*
      \_______/
       ________  <-- gap
      /        \ <-- corner diamond
      |        |
      |    O   | <-- ceter hole
      |        |
      \________/
      */

      const Vector::Triple xyz = (posOnSphere - fOrigin).GetCoordinates(fTelCS);
      const double& x = xyz.get<0>();

      // mirror height limit
      if (abs(x) > nSeg*mirrorSegmentHeight) {
        lost = true;
      } else {
        const double& y = xyz.get<1>();
        const double& z = -xyz.get<2>();
        const double phi = fRCurv * atan2(y, z);
        // mirror width limit
        if (abs(phi) > nSeg*mirrorSegmentAngWidth) {
          lost = true;
        } else {
          const double h = abs(x - int(x/mirrorSegmentHeight) * mirrorSegmentHeight);
          // horizontal gaps between mirror segments
          if (h < mirrorGapWidth)
            lost = true;
          else {
            const double w = abs(phi - int(phi/mirrorSegmentAngWidth) * mirrorSegmentAngWidth);
            // vertical gaps between mirror segments
            if (w < angGapWidth)
              lost = true;
            else {
              const double dh = h - 0.5*mirrorSegmentHeight;
              const double dw = w - 0.5*mirrorSegmentHeight;
              const double holeRadius = 10 * mm; // see technical drawing
              // central holes in the segments
              if (dh*dh + dw*dw < holeRadius*holeRadius)
                lost = true;
              else {
                const double xe = 0.5 * mirrorSegmentHeight - abs(dh);
                const double ye = 0.5 * mirrorSegmentHeight - abs(dw);
                const double d = 20 * mm; // see technical drawing
                // corner of the segement
                if (xe + ye < d) 
                  lost = true;
              }
            }
          } 
        }
      } 
    } else { // czech mirrors
      const double mirrorGapWidth =  5*mm;  // educated guess
      const double angGapWidth = 2 * fRCurv * asin(0.5*mirrorGapWidth / fRCurv);  // approx 10 deg at equator
      const double mirrorSegmentHeight = 624.7*mm; //approx
      const double mirrorSegmentHalfWidth = 0.5 * 541*mm; //approx
      const double mirrorSegmentAngHalfWidth = 2 * fRCurv * asin(0.5*mirrorSegmentHalfWidth / fRCurv);
      //const double mirrorSegmentAngHalfWidth = 4.53952*deg;

      const Vector::Triple xyz = (posOnSphere - fOrigin).GetCoordinates(fTelCS);
      const double& x = xyz.get<0>();
      const double& y = xyz.get<1>();
      const double& z = -xyz.get<2>();
      const double phi = fRCurv * atan2(y, z);

      const vector<double> segEle {-3.98*deg, 11.9316*deg, -19.9279*deg, 28.0239*deg};
      const vector<double> segSideHight {((488.15 - 155.15) / 2) * mm, ((467.87 - 143.99) / 2) * mm, 
                                         ((466.56 - 151.11) / 2) * mm, ((472.67 - 133.77) / 2) * mm};

      // mirror height limit
      if (!lost && abs(x) > 3*mirrorSegmentHeight)
        lost = true;
      // mirror width limit
      if (!lost && abs(phi) > 8.3*mirrorSegmentAngHalfWidth)
        lost = true;
      
      // vertical gaps between mirror segments
      if (!lost) {
        const double w = abs(phi - int(phi/mirrorSegmentAngHalfWidth) * mirrorSegmentAngHalfWidth);
        if (w < angGapWidth) {
          const bool odd = bool(int(phi/mirrorSegmentAngHalfWidth)%2);
          const int s = odd ? 1 : -1;
          for (unsigned int i = 0; i < 4; ++i) {
            const double theta = fRCurv * (atan2(-x, z) + s * segEle.at(i));
            if (abs(theta) < segSideHight.at(i))
              lost = true;
          }
        }
      }
 
      // central holes in segemnts
      if (!lost) {
        const bool odd = bool(int(phi/mirrorSegmentAngHalfWidth)%2);
        const int s = odd ? 1 : -1;
        const double holeRadius = 15*mm;
        for (unsigned int i = 0; i < 4; ++i) {
          const double dPhi_halfSegm = 4.53952*deg;
          const double dPhi0 = odd ? dPhi_halfSegm : 0*deg;
          for (unsigned int j = 0; j < 4; ++j) {
            const double thetaA = fRCurv * (atan2(-x, z) - s * segEle.at(i));
            double phiA = fRCurv * (atan2(y, z) + j * 2 * dPhi_halfSegm + dPhi0);
            if (thetaA*thetaA + phiA*phiA < holeRadius*holeRadius)
              lost = true;
            phiA = fRCurv * (atan2(y, z) - j * 2 * dPhi_halfSegm - dPhi0);
            if (thetaA*thetaA + phiA*phiA < holeRadius*holeRadius)
              lost = true;
          }
        }
      }

      // 'trianlgular' gaps between the hexagons
      if (!lost) {
        
      }

    }

    if (lost) {
      if (HARDDEBUG) {
        cout << posOnSphere.GetX(fTelCS) / mm << ','
             << posOnSphere.GetY(fTelCS) / mm << ','
             << posOnSphere.GetZ(fTelCS) / mm << '\n';
      }
      return eLost;
    }
  } // end mirror gaps

  // the deformed mirror is/will not be used -> kick out some day
  if (fMirrorSegmentSigma > 0 || fMirrorRadiusSigma > 0) {
    const MirrorSegment& ms = GetMirrorSegment(posOnSphere);
    const int sphereStatus2 = RTFunctions::Sphere(ms.origin, ms.radius, photonIn, intersections);
    if (sphereStatus2 <= 0) {
      photonOut = photonIn;
      ostringstream err;
      err << "Missed MirrorSegment : status=" << sphereStatus2 << " index=" << ms.index << " !!!";
      ERROR(err);
      return eLost;
    }
    photonOut = intersections.back().first;
    normal = intersections.back().second;
  }

  // CAREFUL!!! mirror reflectivity is NOT 100%!!!
  // study of mirror test documentation showed
  // mean of 88% total reflectivity, 90% of this light in spot
  // At the moment, photons are absorbed but also part of the light goes to the diffuse part
  // CHECK!!!

  const double reflectedOrNot = RandFlat::shoot(&fRandom.GetEngine(), 0, 1);
  const double photonWL = photonOut.GetWavelength();
  if (reflectedOrNot > fReflectivity.Y(photonWL)) {
    return eLost;  // absorbed
  } else {
    // reflect photon on mirror
    RTFunctions::Reflection(photonOut, normal, photonOut);

    if (fMirrorDiffusionTop && fMirrorDiffusionBot) {
      // use light distribution from mirror measurements in Nov 2014 (results from Lenka)
      // only do diffusion if diffusion function is given in xml AND fDoNoHalo==0, aka the TabulatedFunction is not NULL
      const Vector mirrordiffusion =
        RTFunctions::MirrorDiffusion(fRandom, photonOut.GetDirection(), fMirrorDiffusionTop, fMirrorDiffusionBot, verticalPosOnMirror);
      photonOut.SetDirection(mirrordiffusion);
    }

    if (fMirrorAbsorptionTop > 0 || fMirrorAbsorptionBot > 0)
      RTFunctions::Absorption(photonOut, fMirrorAbsorptionTop, fMirrorAbsorptionBot, verticalPosOnMirror, photonOut, fMirrorSize);

    return eCamera;
  }
}


Mirror::MirrorSegment&
Mirror::GetMirrorSegment(const Point& p)
{
  // If there is a per mirror-segment random component
  // First: index mirror-segments
  const double segmentWidth = gMirrorWidth / gNSegments;
  const double x = p.GetX(fTelCS);
  const double y = p.GetY(fTelCS);
  const int ix = x / segmentWidth - (x < 0 ? 1 : 0);
  const int iy = y / segmentWidth - (y < 0 ? 1 : 0);
  const int index = ix * 1000 + iy;  // arbitrary but unique mirror-segment index

  if (!fSegments.count(index)) {  // if this mirror-segment is new -> create it

    const double radius = RandGauss::shoot(&fRandom.GetEngine(), fRCurv, fMirrorRadiusSigma);
    const double centerSegmentX = (0.5 + ix) * segmentWidth;
    const double centerSegmentY = (0.5 + iy) * segmentWidth;
    const double centerSegmentZ = -sqrt(pow(fRCurv, 2) - pow(centerSegmentX, 2) - pow(centerSegmentY, 2));
    const Point centerSegment(centerSegmentX, centerSegmentY, centerSegmentZ, fTelCS);
    const Vector dirSegment(centerSegment - fOrigin);

    const double elevationSegment = atan(centerSegmentX / fRCurv);

    CoordinateSystemPtr aux1CS = CoordinateSystem::Translation(dirSegment, fTelCS);
    CoordinateSystemPtr aux2CS = CoordinateSystem::RotationY(-elevationSegment, aux1CS);
    CoordinateSystemPtr aux3CS = CoordinateSystem::RotationX((centerSegmentY < 0 ? -1 : 1) * (-dirSegment).GetTheta(aux2CS), aux2CS);

    const double phi = RandFlat::shoot(&fRandom.GetEngine(), 0, 360*deg);
    const double theta = acos(RandFlat::shoot(&fRandom.GetEngine(), cos(fMirrorSegmentSigma), 1));
    //const double theta = acos(cos(fMirrorSegmentSigma));

    CoordinateSystemPtr aux4CS = CoordinateSystem::RotationZ(phi, aux3CS);
    CoordinateSystemPtr aux5CS = CoordinateSystem::RotationY(theta, aux4CS);
    CoordinateSystemPtr msCS = CoordinateSystem::RotationZ(-phi, aux5CS);

    //const Vector delta(sinTheta*cos(phi), sinTheta*sin(phi), cosTheta, aux3CS);
    const Vector delta(0,0,1, msCS);
    const Point originPrime(centerSegment + delta * radius);

    //cout << " creating " << index << " " << phi/deg << " " << theta/deg << " " << fMirrorSegmentSigma/deg << " " << cos(fMirrorSegmentSigma) << endl;

    fSegments[index] = MirrorSegment(index, originPrime, radius, msCS);
  }
  return fSegments[index];
}


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

  const int color = 1;
  const int size = 2;

  const double width = gMirrorWidth; // total width
  const int nSeg = gNSegments;
  const int nSubSeg = 3; // sub structure mirror segments... (at least 1!)

  const double segWidth = width / nSeg;

  //cout << " width=" <<width << " segWidth=" << segWidth << endl;

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

    for (int j = 0; j < nSeg; ++j) {

      //const int color = 1 + j + i*nSeg;

      const double x_mir = -gMirrorWidth / 2 + (0.5 + i) * segWidth;
      const double y_mir = -gMirrorWidth / 2 + (0.5 + j) * segWidth;

      //cout << " draw: " << i << j << endl;
      const MirrorSegment& ms = GetMirrorSegment(Point(x_mir, y_mir, 0, fTelCS));
      const double r = ms.radius;

      const double d = segWidth / nSubSeg;

      /*cout << " draw: " << i << " " << j << " " << x_mir << " " << y_mir << " "
           << " index=" << ms.index << " r=" << ms.radius << " "
           << ms.origin.GetX(fTelCS)<< " "
           << ms.origin.GetY(fTelCS)<< " "
           << ms.origin.GetZ(fTelCS)<< " "
           <<  endl;*/

      for (int k = 0; k < nSubSeg; ++k) {

        for (int l = 0; l < nSubSeg; ++l) {

          const double x1 = -segWidth / 2 + k * d;
          const double y1 = -segWidth / 2 + l * d;
          const double z1 = r - sqrt(pow(r, 2) - pow(x1, 2) - pow(y1, 2));
          const Point p1(x1, y1, z1, ms.cs);

          //cout << " z1=" << z1 << " " << p1.GetZ(ms.cs) << " " << p1.GetZ(fTelCS) << endl;

          const double x2 = -segWidth / 2 + (k + 1) * d;
          const double y2 = -segWidth / 2 + l * d;
          const double z2 = r - sqrt(pow(r, 2) - pow(x2, 2) - pow(y2, 2));
          const Point p2(x2, y2, z2, ms.cs);

          const double x3 = -segWidth / 2 + k * d;
          const double y3 = -segWidth / 2 + (l + 1) * d;
          const double z3 = r - sqrt(pow(r, 2) - pow(x3, 2) - pow(y3, 2));
          const Point p3(x3, y3, z3, ms.cs);

          const double x4 = -segWidth / 2 + (k + 1) * d;
          const double y4 = -segWidth / 2 + (l + 1) * d;
          const double z4 = r - sqrt(pow(r, 2) - pow(x4, 2) - pow(y4, 2));
          const Point p4(x4, y4, z4, ms.cs);

          TPolyLine3D* const l1 = new TPolyLine3D(5);
          l1->SetLineColor(color);
          l1->SetLineWidth(size);
          l1->SetPoint(0, p1.GetX(fTelCS), p1.GetY(fTelCS), p1.GetZ(fTelCS));
          l1->SetPoint(1, p2.GetX(fTelCS), p2.GetY(fTelCS), p2.GetZ(fTelCS));
          l1->SetPoint(2, p4.GetX(fTelCS), p4.GetY(fTelCS), p4.GetZ(fTelCS));
          l1->SetPoint(3, p3.GetX(fTelCS), p3.GetY(fTelCS), p3.GetZ(fTelCS));
          l1->SetPoint(4, p1.GetX(fTelCS), p1.GetY(fTelCS), p1.GetZ(fTelCS));

          l1->Draw();
          objs->AddLast(l1);
        }
      }
    }
  }
  return objs;
}