TEveEllipsoidGL.cc

Go to the documentation of this file.

#include "Fireworks/Vertices/interface/TEveEllipsoidGL.h"

#include "TEveProjections.h"  // AMT missing getter for projection center / beam-spot

#include "Fireworks/Vertices/interface/TEveEllipsoidGL.h"
#include "Fireworks/Vertices/interface/TEveEllipsoid.h"
#include "TEveProjectionManager.h"

#include "TMath.h"

#include "TGLRnrCtx.h"

#include "TMatrixDEigen.h"
#include "TMatrixDSym.h"

#include "TDecompSVD.h"
#include "TGLIncludes.h"

//==============================================================================
// TEveEllipsoidGL
//==============================================================================

//______________________________________________________________________________
// OpenGL renderer class for TEveEllipsoid.
//

//______________________________________________________________________________
TEveEllipsoidGL::TEveEllipsoidGL() : TGLObject(), fE(nullptr) {
  // Constructor.

  // fDLCache = kFALSE; // Disable display list.
}

//______________________________________________________________________________
Bool_t TEveEllipsoidGL::SetModel(TObject* obj, const Option_t* /*opt*/) {
  // Set model object.

  fE = SetModelDynCast<TEveEllipsoid>(obj);
  return kTRUE;
}

//______________________________________________________________________________
void TEveEllipsoidGL::SetBBox() {
  // Set bounding box.
  ((TEveEllipsoid*)fExternalObj)->ComputeBBox();
  SetAxisAlignedBBox(((TEveEllipsoid*)fExternalObj)->AssertBBox());
}

namespace {
  GLUquadric* quad = nullptr;  // !!!! AMT check why TGLQuadric crashes on mac
}

//______________________________________________________________________________
void TEveEllipsoidGL::DirectDraw(TGLRnrCtx& /*rnrCtx*/) const {
  // Render with OpenGL.

  // printf("TEveEllipsoidGL::DirectDraw LOD %s\n", fE->GetName());

  glPushAttrib(GL_ENABLE_BIT | GL_POLYGON_BIT | GL_LIGHTING_BIT);
  glEnable(GL_NORMALIZE);
  if (!quad)
    quad = gluNewQuadric();

  glPushMatrix();

  TMatrixDSym xxx(3);
  for (int i = 0; i < 3; i++)
    for (int j = 0; j < 3; j++) {
      xxx(i, j) = fE->RefEMtx()(i + 1, j + 1);
    }
  TMatrixDEigen eig(xxx);

  // rewrite for multmatrix ....
  TEveTrans x;
  for (int i = 0; i < 3; i++)
    for (int j = 0; j < 3; j++) {
      x(i + 1, j + 1) = eig.GetEigenVectors()(i, j);
    }

  TVector3 a = x.GetBaseVec(1);
  TVector3 c = a.Cross(x.GetBaseVec(2));
  x.SetBaseVec(3, c);

  glTranslatef(fE->RefPos()[0], fE->RefPos()[1], fE->RefPos()[2]);
  glMultMatrixd(x.Array());
  glScalef(fE->RefExtent3D()[0], fE->RefExtent3D()[1], fE->RefExtent3D()[2]);
  gluSphere(quad, 1., 30, 30);

  glPopMatrix();
  glPopAttrib();

  // gluDeleteQuadric(quad);
}

//==============================================================================
// TEveEllipsoidProjectedGL
//==============================================================================

//______________________________________________________________________________
// OpenGL renderer class for TEveEllipsoidProjected.
//

//______________________________________________________________________________
TEveEllipsoidProjectedGL::TEveEllipsoidProjectedGL() : fM(nullptr) {
  // Constructor.

  // fDLCache = kFALSE; // Disable display list.
  fMultiColor = kTRUE;
}

//______________________________________________________________________________
Bool_t TEveEllipsoidProjectedGL::SetModel(TObject* obj, const Option_t* /*opt*/) {
  // Set model object.

  fM = SetModelDynCast<TEveEllipsoidProjected>(obj);
  fE = dynamic_cast<TEveEllipsoid*>(fM->GetProjectable());
  return fE != nullptr;
}

//______________________________________________________________________________
void TEveEllipsoidProjectedGL::SetBBox() {
  // Set bounding box.

  SetAxisAlignedBBox(((TEveEllipsoidProjected*)fExternalObj)->AssertBBox());
}

//______________________________________________________________________________
void TEveEllipsoidProjectedGL::DirectDraw(TGLRnrCtx& rnrCtx) const {
  // Render with OpenGL.

  TEveProjection* proj = fM->GetManager()->GetProjection();

  glPushAttrib(GL_ENABLE_BIT | GL_POLYGON_BIT | GL_LINE_BIT | GL_POINT_BIT);
  glDisable(GL_LIGHTING);
  glDisable(GL_CULL_FACE);

  glPushMatrix();
  if (proj->GetType() == TEveProjection::kPT_RPhi)
    DrawRhoPhi();
  else
    DrawRhoZ();

  glPopMatrix();
  glPopAttrib();
}

//______________________________________________________________________________
void TEveEllipsoidProjectedGL::drawArch(
    float phiStart, float phiEnd, float phiStep, TEveVector& v0, TEveVector& v1, TEveVector& v2) const {
  TEveProjection* proj = fM->GetManager()->GetProjection();
  float phi = phiStart;
  while (phi < phiEnd) {
    TEveVector v = v0 + v1 * ((float)cos(phi)) + v2 * ((float)sin(phi));
    proj->ProjectVector(v, fM->fDepth);
    glVertex3fv(v.Arr());

    phi += phiStep;
  }
  TEveVector v = v0 + v1 * ((float)cos(phiEnd)) + v2 * ((float)sin(phiEnd));
  proj->ProjectVector(v, fM->fDepth);
  glVertex3fv(v.Arr());
}

//-------------------------------------------------------------------------------
void TEveEllipsoidProjectedGL::DrawRhoPhi() const {
  // printf("TEveEllipsoidProjectedGL::DirectDraw [%s ]\n", fE->GetName() );

  TMatrixDSym xxx(3);
  for (int i = 0; i < 2; i++)
    for (int j = 0; j < 2; j++) {
      xxx(i, j) = fE->RefEMtx()(i + 1, j + 1);
    }

  TMatrixDEigen eig(xxx);
  TVectorD xxxEig(eig.GetEigenValuesRe());

  // Projection supports only floats  :(
  TEveVector v0(fE->RefPos()[0], fE->RefPos()[1], 0);
  TEveVector v1(eig.GetEigenVectors()(0, 0), eig.GetEigenVectors()(0, 1), 0);
  v1 *= fE->fEScale * sqrt(TMath::Abs(xxxEig[0]));
  TEveVector v2(eig.GetEigenVectors()(1, 0), eig.GetEigenVectors()(1, 1), 0);
  v2 *= fE->fEScale * sqrt(TMath::Abs(xxxEig[1]));

  TEveProjection* proj = fM->GetManager()->GetProjection();

  // fill
  glBegin(GL_POLYGON);
  drawArch(0, TMath::TwoPi(), TMath::TwoPi() / 20, v0, v1, v2);
  glEnd();

  // frame
  TGLUtil::LineWidth(fE->fLineWidth);
  TGLUtil::Color(fE->fLineColor);

  glBegin(GL_LINE_LOOP);
  drawArch(0, TMath::TwoPi(), TMath::TwoPi() / 20, v0, v1, v2);
  glEnd();

  glBegin(GL_LINES);
  {
    // glColor3f(1, 0, 0);
    TEveVector p1 = v0 - v1;
    TEveVector p2 = v0 + v1;
    proj->ProjectVector(p1, fM->fDepth);
    proj->ProjectVector(p2, fM->fDepth);
    glVertex3fv(p1.Arr());
    glVertex3fv(p2.Arr());
  }
  {
    //  glColor3f(0, 1, 0);
    TEveVector p1 = v0 - v2;
    TEveVector p2 = v0 + v2;
    proj->ProjectVector(p1, fM->fDepth);
    proj->ProjectVector(p2, fM->fDepth);
    glVertex3fv(p1.Arr());
    glVertex3fv(p2.Arr());
  }
  glEnd();
}

//--------------------------------------------------------------------
void TEveEllipsoidProjectedGL::DrawRhoZ() const {
  // printf("TEveEllipsoidProjectedGL::DirectDraw [%s ]\n", fE->GetTitle() );

  TEveVector v0(fE->RefPos()[0], fE->RefPos()[1], fE->RefPos()[2]);

  TMatrixDSym xxx(3);
  xxx(0, 0) = 1;
  for (int i = 1; i < 3; i++)
    for (int j = 1; j < 3; j++) {
      xxx(i, j) = fE->RefEMtx()(i + 1, j + 1);
    }
  TMatrixDEigen eig(xxx);
  TVectorD xxxEig(eig.GetEigenValuesRe());

  TEveVector v1(0, eig.GetEigenVectors()(1, 2), eig.GetEigenVectors()(2, 2));
  v1 *= fE->fEScale * sqrt(TMath::Abs(xxxEig[2]));

  TEveVector v2(0, eig.GetEigenVectors()(1, 1), eig.GetEigenVectors()(2, 1));
  v2 *= fE->fEScale * sqrt(TMath::Abs(xxxEig[1]));
  if (v1[1] * v2[2] > v1[2] * v2[1])
    v2 *= -1;

  TEveProjection* proj = fM->GetManager()->GetProjection();

  // ellipse intersection with projection center
  bool splitted = false;
  int N = 20;
  double phiStep = TMath::TwoPi() / N;

  // projection center can be moved in beam-spot
  float bs = 0;
  if (proj->GetDisplaceOrigin())
    bs = proj->RefCenter()[1];

  float da = v2[1] * v2[1] + v1[1] * v1[1];
  float db = 2 * v1[1] * (v0[1] - bs);
  float dc = (v0[1] - bs) * (v0[1] - bs) - v2[1] * v2[1];
  float disc = (db * db - 4 * da * dc);

  if (disc > 0) {
    disc = sqrt(disc);
    float cosS1 = (-db + disc) / (2 * da);
    float cosS2 = (-db - disc) / (2 * da);
    if (TMath::Abs(cosS1) < 1) {
      splitted = true;
      // printf("splitted \n");

      double phi1 = acos(cosS1);
      double phi2 = acos(cosS2);
      TEveVector ps1 = v0 + v1 * ((float)cos(phi1)) + v2 * ((float)sin(phi1));
      TEveVector ps2 = v0 + v1 * ((float)cos(phi2)) + v2 * ((float)sin(phi2));

      // acos has values [0, Pi] , check the symetry over x axis (mirroring)
      if (TMath::Abs(ps1[1] - bs) > 1e-5)
        phi1 = TMath::TwoPi() - phi1;

      if (TMath::Abs(ps2[1] - bs) > 1e-5)
        phi2 = TMath::TwoPi() - phi2;

      int N = 20;
      double phiStep = TMath::TwoPi() / N;
      double phiOffset = phiStep * 0.1;
      double phiMin = TMath::Min(phi1, phi2);
      double phiMax = TMath::Max(phi1, phi2);
      // printf(" %f %f \n",phi1*TMath::RadToDeg(), phi2*TMath::RadToDeg() );

      // fill
      // upper clothing
      glBegin(GL_POLYGON);
      drawArch(phiMin + phiOffset, phiMax - phiOffset, phiStep, v0, v1, v2);
      glEnd();
      // bottom clothing
      glBegin(GL_POLYGON);
      drawArch(phiMax + phiOffset, phiMax + TMath::TwoPi() - (phiMax - phiMin) - phiOffset, phiStep, v0, v1, v2);
      glEnd();

      // frame
      TGLUtil::LineWidth(fE->fLineWidth);
      TGLUtil::Color(fE->fLineColor);
      // upper clothing
      glBegin(GL_LINE_LOOP);
      drawArch(phiMin + phiOffset, phiMax - phiOffset, phiStep, v0, v1, v2);
      glEnd();
      // bottom clothing
      glBegin(GL_LINE_LOOP);
      drawArch(phiMax + phiOffset, phiMax + TMath::TwoPi() - (phiMax - phiMin) - phiOffset, phiStep, v0, v1, v2);
      glEnd();
    }
  }

  if (!splitted) {
    glBegin(GL_POLYGON);
    drawArch(0, TMath::TwoPi(), phiStep, v0, v1, v2);
    glEnd();
    TGLUtil::LineWidth(fE->fLineWidth);
    TGLUtil::Color(fE->fLineColor);
    glBegin(GL_LINE_LOOP);
    drawArch(0, TMath::TwoPi(), phiStep, v0, v1, v2);
    glEnd();
  }

  drawRhoZAxis(v0, v2);
  drawRhoZAxis(v0, v1);
}

//______________________________________________________________________________
void TEveEllipsoidProjectedGL::drawRhoZAxis(TEveVector& v0, TEveVector& v2) const {
  glBegin(GL_LINES);
  TEveProjection* proj = fM->GetManager()->GetProjection();

  float bs = 0;
  if (proj->GetDisplaceOrigin())
    bs = proj->RefCenter()[1];

  float off = (v2[1] > v0[1]) ? 0.01 : -0.01;
  TEveVector alu = v0 + v2;
  proj->ProjectVector(alu, fM->fDepth);
  glVertex3fv(alu.Arr());

  if (TMath::Abs(v0[1] / v2[1]) < 1) {
    alu = v0 - ((float)((1 - off) * (v0[1] - bs) / v2[1])) * v2;
    proj->ProjectVector(alu, fM->fDepth);
    glVertex3fv(alu.Arr());

    //============================

    alu = v0 - ((float)((1 + off) * (v0[1] - bs) / v2[1])) * v2;
    proj->ProjectVector(alu, fM->fDepth);
    glVertex3fv(alu.Arr());
  }

  alu = v0 - v2;
  proj->ProjectVector(alu, fM->fDepth);
  glVertex3fv(alu.Arr());

  glEnd();
}