Vertices/src/TEveEllipsoidGL.cc

0001 #include "Fireworks/Vertices/interface/TEveEllipsoidGL.h"
0002
0003 #include "TEveProjections.h"  // AMT missing getter for projection center / beam-spot
0004
0005 #include "Fireworks/Vertices/interface/TEveEllipsoidGL.h"
0006 #include "Fireworks/Vertices/interface/TEveEllipsoid.h"
0007 #include "TEveProjectionManager.h"
0008
0009 #include "TMath.h"
0010
0011 #include "TGLRnrCtx.h"
0012
0013 #include "TMatrixDEigen.h"
0014 #include "TMatrixDSym.h"
0015
0016 #include "TDecompSVD.h"
0017 #include "TGLIncludes.h"
0018
0019 //==============================================================================
0020 // TEveEllipsoidGL
0021 //==============================================================================
0022
0023 //______________________________________________________________________________
0024 // OpenGL renderer class for TEveEllipsoid.
0025 //
0026
0027 //______________________________________________________________________________
0028 TEveEllipsoidGL::TEveEllipsoidGL() : TGLObject(), fE(nullptr) {
0029   // Constructor.
0030
0031   // fDLCache = kFALSE; // Disable display list.
0032 }
0033
0034 //______________________________________________________________________________
0035 Bool_t TEveEllipsoidGL::SetModel(TObject* obj, const Option_t* /*opt*/) {
0036   // Set model object.
0037
0038   fE = SetModelDynCast<TEveEllipsoid>(obj);
0039   return kTRUE;
0040 }
0041
0042 //______________________________________________________________________________
0043 void TEveEllipsoidGL::SetBBox() {
0044   // Set bounding box.
0045   ((TEveEllipsoid*)fExternalObj)->ComputeBBox();
0046   SetAxisAlignedBBox(((TEveEllipsoid*)fExternalObj)->AssertBBox());
0047 }
0048
0049 namespace {
0050   GLUquadric* quad = nullptr;  // !!!! AMT check why TGLQuadric crashes on mac
0051 }
0052
0053 //______________________________________________________________________________
0054 void TEveEllipsoidGL::DirectDraw(TGLRnrCtx& /*rnrCtx*/) const {
0055   // Render with OpenGL.
0056
0057   // printf("TEveEllipsoidGL::DirectDraw LOD %s\n", fE->GetName());
0058
0059   glPushAttrib(GL_ENABLE_BIT | GL_POLYGON_BIT | GL_LIGHTING_BIT);
0060   glEnable(GL_NORMALIZE);
0061   if (!quad)
0062     quad = gluNewQuadric();
0063
0064   glPushMatrix();
0065
0066   TMatrixDSym xxx(3);
0067   for (int i = 0; i < 3; i++)
0068     for (int j = 0; j < 3; j++) {
0069       xxx(i, j) = fE->RefEMtx()(i + 1, j + 1);
0070     }
0071   TMatrixDEigen eig(xxx);
0072
0073   // rewrite for multmatrix ....
0074   TEveTrans x;
0075   for (int i = 0; i < 3; i++)
0076     for (int j = 0; j < 3; j++) {
0077       x(i + 1, j + 1) = eig.GetEigenVectors()(i, j);
0078     }
0079
0080   TVector3 a = x.GetBaseVec(1);
0081   TVector3 c = a.Cross(x.GetBaseVec(2));
0082   x.SetBaseVec(3, c);
0083
0084   glTranslatef(fE->RefPos()[0], fE->RefPos()[1], fE->RefPos()[2]);
0085   glMultMatrixd(x.Array());
0086   glScalef(fE->RefExtent3D()[0], fE->RefExtent3D()[1], fE->RefExtent3D()[2]);
0087   gluSphere(quad, 1., 30, 30);
0088
0089   glPopMatrix();
0090   glPopAttrib();
0091
0092   // gluDeleteQuadric(quad);
0093 }
0094
0095 //==============================================================================
0096 // TEveEllipsoidProjectedGL
0097 //==============================================================================
0098
0099 //______________________________________________________________________________
0100 // OpenGL renderer class for TEveEllipsoidProjected.
0101 //
0102
0103 //______________________________________________________________________________
0104 TEveEllipsoidProjectedGL::TEveEllipsoidProjectedGL() : fM(nullptr) {
0105   // Constructor.
0106
0107   // fDLCache = kFALSE; // Disable display list.
0108   fMultiColor = kTRUE;
0109 }
0110
0111 //______________________________________________________________________________
0112 Bool_t TEveEllipsoidProjectedGL::SetModel(TObject* obj, const Option_t* /*opt*/) {
0113   // Set model object.
0114
0115   fM = SetModelDynCast<TEveEllipsoidProjected>(obj);
0116   fE = dynamic_cast<TEveEllipsoid*>(fM->GetProjectable());
0117   return fE != nullptr;
0118 }
0119
0120 //______________________________________________________________________________
0121 void TEveEllipsoidProjectedGL::SetBBox() {
0122   // Set bounding box.
0123
0124   SetAxisAlignedBBox(((TEveEllipsoidProjected*)fExternalObj)->AssertBBox());
0125 }
0126
0127 //______________________________________________________________________________
0128 void TEveEllipsoidProjectedGL::DirectDraw(TGLRnrCtx& rnrCtx) const {
0129   // Render with OpenGL.
0130
0131   TEveProjection* proj = fM->GetManager()->GetProjection();
0132
0133   glPushAttrib(GL_ENABLE_BIT | GL_POLYGON_BIT | GL_LINE_BIT | GL_POINT_BIT);
0134   glDisable(GL_LIGHTING);
0135   glDisable(GL_CULL_FACE);
0136
0137   glPushMatrix();
0138   if (proj->GetType() == TEveProjection::kPT_RPhi)
0139     DrawRhoPhi();
0140   else
0141     DrawRhoZ();
0142
0143   glPopMatrix();
0144   glPopAttrib();
0145 }
0146
0147 //______________________________________________________________________________
0148 void TEveEllipsoidProjectedGL::drawArch(
0149     float phiStart, float phiEnd, float phiStep, TEveVector& v0, TEveVector& v1, TEveVector& v2) const {
0150   TEveProjection* proj = fM->GetManager()->GetProjection();
0151   float phi = phiStart;
0152   while (phi < phiEnd) {
0153     TEveVector v = v0 + v1 * ((float)cos(phi)) + v2 * ((float)sin(phi));
0154     proj->ProjectVector(v, fM->fDepth);
0155     glVertex3fv(v.Arr());
0156
0157     phi += phiStep;
0158   }
0159   TEveVector v = v0 + v1 * ((float)cos(phiEnd)) + v2 * ((float)sin(phiEnd));
0160   proj->ProjectVector(v, fM->fDepth);
0161   glVertex3fv(v.Arr());
0162 }
0163
0164 //-------------------------------------------------------------------------------
0165 void TEveEllipsoidProjectedGL::DrawRhoPhi() const {
0166   // printf("TEveEllipsoidProjectedGL::DirectDraw [%s ]\n", fE->GetName() );
0167
0168   TMatrixDSym xxx(3);
0169   for (int i = 0; i < 2; i++)
0170     for (int j = 0; j < 2; j++) {
0171       xxx(i, j) = fE->RefEMtx()(i + 1, j + 1);
0172     }
0173
0174   TMatrixDEigen eig(xxx);
0175   TVectorD xxxEig(eig.GetEigenValuesRe());
0176
0177   // Projection supports only floats  :(
0178   TEveVector v0(fE->RefPos()[0], fE->RefPos()[1], 0);
0179   TEveVector v1(eig.GetEigenVectors()(0, 0), eig.GetEigenVectors()(0, 1), 0);
0180   v1 *= fE->fEScale * sqrt(TMath::Abs(xxxEig[0]));
0181   TEveVector v2(eig.GetEigenVectors()(1, 0), eig.GetEigenVectors()(1, 1), 0);
0182   v2 *= fE->fEScale * sqrt(TMath::Abs(xxxEig[1]));
0183
0184   TEveProjection* proj = fM->GetManager()->GetProjection();
0185
0186   // fill
0187   glBegin(GL_POLYGON);
0188   drawArch(0, TMath::TwoPi(), TMath::TwoPi() / 20, v0, v1, v2);
0189   glEnd();
0190
0191   // frame
0192   TGLUtil::LineWidth(fE->fLineWidth);
0193   TGLUtil::Color(fE->fLineColor);
0194
0195   glBegin(GL_LINE_LOOP);
0196   drawArch(0, TMath::TwoPi(), TMath::TwoPi() / 20, v0, v1, v2);
0197   glEnd();
0198
0199   glBegin(GL_LINES);
0200   {
0201     // glColor3f(1, 0, 0);
0202     TEveVector p1 = v0 - v1;
0203     TEveVector p2 = v0 + v1;
0204     proj->ProjectVector(p1, fM->fDepth);
0205     proj->ProjectVector(p2, fM->fDepth);
0206     glVertex3fv(p1.Arr());
0207     glVertex3fv(p2.Arr());
0208   }
0209   {
0210     //  glColor3f(0, 1, 0);
0211     TEveVector p1 = v0 - v2;
0212     TEveVector p2 = v0 + v2;
0213     proj->ProjectVector(p1, fM->fDepth);
0214     proj->ProjectVector(p2, fM->fDepth);
0215     glVertex3fv(p1.Arr());
0216     glVertex3fv(p2.Arr());
0217   }
0218   glEnd();
0219 }
0220
0221 //--------------------------------------------------------------------
0222 void TEveEllipsoidProjectedGL::DrawRhoZ() const {
0223   // printf("TEveEllipsoidProjectedGL::DirectDraw [%s ]\n", fE->GetTitle() );
0224
0225   TEveVector v0(fE->RefPos()[0], fE->RefPos()[1], fE->RefPos()[2]);
0226
0227   TMatrixDSym xxx(3);
0228   xxx(0, 0) = 1;
0229   for (int i = 1; i < 3; i++)
0230     for (int j = 1; j < 3; j++) {
0231       xxx(i, j) = fE->RefEMtx()(i + 1, j + 1);
0232     }
0233   TMatrixDEigen eig(xxx);
0234   TVectorD xxxEig(eig.GetEigenValuesRe());
0235
0236   TEveVector v1(0, eig.GetEigenVectors()(1, 2), eig.GetEigenVectors()(2, 2));
0237   v1 *= fE->fEScale * sqrt(TMath::Abs(xxxEig[2]));
0238
0239   TEveVector v2(0, eig.GetEigenVectors()(1, 1), eig.GetEigenVectors()(2, 1));
0240   v2 *= fE->fEScale * sqrt(TMath::Abs(xxxEig[1]));
0241   if (v1[1] * v2[2] > v1[2] * v2[1])
0242     v2 *= -1;
0243
0244   TEveProjection* proj = fM->GetManager()->GetProjection();
0245
0246   // ellipse intersection with projection center
0247   bool splitted = false;
0248   int N = 20;
0249   double phiStep = TMath::TwoPi() / N;
0250
0251   // projection center can be moved in beam-spot
0252   float bs = 0;
0253   if (proj->GetDisplaceOrigin())
0254     bs = proj->RefCenter()[1];
0255
0256   float da = v2[1] * v2[1] + v1[1] * v1[1];
0257   float db = 2 * v1[1] * (v0[1] - bs);
0258   float dc = (v0[1] - bs) * (v0[1] - bs) - v2[1] * v2[1];
0259   float disc = (db * db - 4 * da * dc);
0260
0261   if (disc > 0) {
0262     disc = sqrt(disc);
0263     float cosS1 = (-db + disc) / (2 * da);
0264     float cosS2 = (-db - disc) / (2 * da);
0265     if (TMath::Abs(cosS1) < 1) {
0266       splitted = true;
0267       // printf("splitted \n");
0268
0269       double phi1 = acos(cosS1);
0270       double phi2 = acos(cosS2);
0271       TEveVector ps1 = v0 + v1 * ((float)cos(phi1)) + v2 * ((float)sin(phi1));
0272       TEveVector ps2 = v0 + v1 * ((float)cos(phi2)) + v2 * ((float)sin(phi2));
0273
0274       // acos has values [0, Pi] , check the symetry over x axis (mirroring)
0275       if (TMath::Abs(ps1[1] - bs) > 1e-5)
0276         phi1 = TMath::TwoPi() - phi1;
0277
0278       if (TMath::Abs(ps2[1] - bs) > 1e-5)
0279         phi2 = TMath::TwoPi() - phi2;
0280
0281       int N = 20;
0282       double phiStep = TMath::TwoPi() / N;
0283       double phiOffset = phiStep * 0.1;
0284       double phiMin = TMath::Min(phi1, phi2);
0285       double phiMax = TMath::Max(phi1, phi2);
0286       // printf(" %f %f \n",phi1*TMath::RadToDeg(), phi2*TMath::RadToDeg() );
0287
0288       // fill
0289       // upper clothing
0290       glBegin(GL_POLYGON);
0291       drawArch(phiMin + phiOffset, phiMax - phiOffset, phiStep, v0, v1, v2);
0292       glEnd();
0293       // bottom clothing
0294       glBegin(GL_POLYGON);
0295       drawArch(phiMax + phiOffset, phiMax + TMath::TwoPi() - (phiMax - phiMin) - phiOffset, phiStep, v0, v1, v2);
0296       glEnd();
0297
0298       // frame
0299       TGLUtil::LineWidth(fE->fLineWidth);
0300       TGLUtil::Color(fE->fLineColor);
0301       // upper clothing
0302       glBegin(GL_LINE_LOOP);
0303       drawArch(phiMin + phiOffset, phiMax - phiOffset, phiStep, v0, v1, v2);
0304       glEnd();
0305       // bottom clothing
0306       glBegin(GL_LINE_LOOP);
0307       drawArch(phiMax + phiOffset, phiMax + TMath::TwoPi() - (phiMax - phiMin) - phiOffset, phiStep, v0, v1, v2);
0308       glEnd();
0309     }
0310   }
0311
0312   if (!splitted) {
0313     glBegin(GL_POLYGON);
0314     drawArch(0, TMath::TwoPi(), phiStep, v0, v1, v2);
0315     glEnd();
0316     TGLUtil::LineWidth(fE->fLineWidth);
0317     TGLUtil::Color(fE->fLineColor);
0318     glBegin(GL_LINE_LOOP);
0319     drawArch(0, TMath::TwoPi(), phiStep, v0, v1, v2);
0320     glEnd();
0321   }
0322
0323   drawRhoZAxis(v0, v2);
0324   drawRhoZAxis(v0, v1);
0325 }
0326
0327 //______________________________________________________________________________
0328 void TEveEllipsoidProjectedGL::drawRhoZAxis(TEveVector& v0, TEveVector& v2) const {
0329   glBegin(GL_LINES);
0330   TEveProjection* proj = fM->GetManager()->GetProjection();
0331
0332   float bs = 0;
0333   if (proj->GetDisplaceOrigin())
0334     bs = proj->RefCenter()[1];
0335
0336   float off = (v2[1] > v0[1]) ? 0.01 : -0.01;
0337   TEveVector alu = v0 + v2;
0338   proj->ProjectVector(alu, fM->fDepth);
0339   glVertex3fv(alu.Arr());
0340
0341   if (TMath::Abs(v0[1] / v2[1]) < 1) {
0342     alu = v0 - ((float)((1 - off) * (v0[1] - bs) / v2[1])) * v2;
0343     proj->ProjectVector(alu, fM->fDepth);
0344     glVertex3fv(alu.Arr());
0345
0346     //============================
0347
0348     alu = v0 - ((float)((1 + off) * (v0[1] - bs) / v2[1])) * v2;
0349     proj->ProjectVector(alu, fM->fDepth);
0350     glVertex3fv(alu.Arr());
0351   }
0352
0353   alu = v0 - v2;
0354   proj->ProjectVector(alu, fM->fDepth);
0355   glVertex3fv(alu.Arr());
0356
0357   glEnd();
0358 }