ZetaPixel.cc

Go to the documentation of this file.

#include "ZetaPixel.h"
#include <utl/MathConstants.h>
#include <utl/AugerUnits.h>
#include <vector>
#include <cmath>
#include <algorithm>
#include <utl/TabulatedFunction.h>

using namespace std;
using namespace utl;
using namespace FdEnergyDepositFinderKG;


//---------------------------------------------------
ZetaPixel::ZetaPixel(const unsigned int pixelId,
                     const double xPos, const double yPos,
                     const double sideLength, const double tiltAngle) :
  fLightFraction(-1),
  fPixelId(pixelId),
  fCenterX(xPos),
  fCenterY(yPos),
  fSideLength(sideLength),
  fTiltAngle(tiltAngle),
  fCosTiltAngle(cos(tiltAngle)),
  fSinTiltAngle(sin(tiltAngle))
{
  CalculateCoordinates();
}


//---------------------------------------------------
double
ZetaPixel::CalculateInsideArcLength(const double circleX,
                                    const double circleY,
                                    const double circleRadius)
  const
{
  vector<double> phiVec = CircleIntersections(circleX, circleY, circleRadius);

  // if there is no intercept it is either all in or all out
  if (phiVec.empty())
    return IsInside(circleX, circleY+circleRadius) ? 1 : 0;

  phiVec.push_back(phiVec.front() + 2*kPi);

  double arcLength = 0.;
  for (unsigned int i = 0, n = phiVec.size()-1; i < n; ++i) {
    const double phi1 = phiVec[i];
    const double phi2 = phiVec[i+1];
    const double meanPhi = 0.5*(phi1 + phi2);
    // check for inside or outside at middle of arc
    if (IsInside(circleX + circleRadius*cos(meanPhi),
                 circleY + circleRadius*sin(meanPhi)))
      arcLength += phi2 - phi1;
  }
  return 0.5*arcLength/kPi;
}


//---------------------------------------------------
vector<double>
ZetaPixel::CircleIntersections(const double circleX,
                               const double circleY,
                               const double circleRadius)
  const
{
  vector<double> intersectionPhi;

  const double r = circleRadius;
  const double r2 = r*r;

  // loop over all sides, calculate circle-line intersections
  // (see mathworld.wolfram.com/Circle-LineIntersection.html)
  // and pick the ones that are on the hexagon

  for (unsigned int iSide = 0; iSide < 6; ++iSide) {

    const double x1 = fXCoordinates[iSide] - circleX;
    const double y1 = fYCoordinates[iSide] - circleY;
    const double x2 = fXCoordinates[iSide+1] - circleX;
    const double y2 = fYCoordinates[iSide+1] - circleY;
    const double dx = x2 - x1;
    const double dy = y2 - y1;
    const double dr2 = dx*dx + dy*dy;
    const double D = x1*y2 - x2*y1;
    const double D2 = D*D;
    const double Delta = r2*dr2 - D2;
    if (Delta > 0) {
      const double signdy = (dy < 0 ? -1 : 1);
      const double sqrtDelta = sqrt(Delta);
      const double tx1 = D*dy;
      const double tx2 = signdy*dx*sqrtDelta;
      const double ty1 = -D*dx;
      const double ty2 = fabs(dy)*sqrtDelta;

      const double xInter1 = (tx1 + tx2)/dr2 + circleX;
      const double xInter2 = (tx1 - tx2)/dr2 + circleX;
      const double yInter1 = (ty1 + ty2)/dr2 + circleY;
      const double yInter2 = (ty1 - ty2)/dr2 + circleY;
      if (IsInside(xInter1, yInter1)) {
        double phi = atan2(yInter1-circleY, xInter1-circleX);
        if (phi < 0)
          phi += 2*kPi;
        intersectionPhi.push_back(phi);
      }
      if (IsInside(xInter2, yInter2)) {
        double phi = atan2(yInter2-circleY, xInter2-circleX);
        if (phi < 0)
          phi += 2*kPi;
        intersectionPhi.push_back(phi);
      }
    }
  }

  sort(intersectionPhi.begin(), intersectionPhi.end());

  return intersectionPhi;
}


//---------------------------------------------------
bool
ZetaPixel::IsInside(const double x, const double y)
  const
{
  const double safety = 1e-10;
  const double cos30 = 0.5*kSqrt3;
  const double tan30 = 1/kSqrt3;

  const double r = (1+safety)*fSideLength;
  const double xx = x-fCenterX;
  const double yy = y-fCenterY;
  const double xxRot = fabs(fCosTiltAngle*xx +fSinTiltAngle*yy);
  const double yyRot = fabs(-fSinTiltAngle*xx+fCosTiltAngle*yy);
  const double xMax = cos30*r;

  // inside hexagon check
  return (xxRot <= xMax && yyRot <= r-tan30*xxRot);
}


//---------------------------------------------------
void
ZetaPixel::CalculateCoordinates()
{
  double phi = 30.*degree + fTiltAngle;
  for (unsigned int i = 0; i < 7; ++i) {
    const double x = cos(phi)*fSideLength + fCenterX;
    const double y = sin(phi)*fSideLength + fCenterY;
    fXCoordinates.push_back(x);
    fYCoordinates.push_back(y);
    phi += 60*degree;
  }
}