// @(#)root/g3d:$Id: TRotMatrix.cxx 20882 2007-11-19 11:31:26Z rdm $
// Author: Rene Brun   14/09/95

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#include "TGeometry.h"
#include "TRotMatrix.h"
#include "TClass.h"
#include "TMath.h"

ClassImp(TRotMatrix)


//______________________________________________________________________________
//
// Manages a detector rotation matrix. See class TGeometry.
//


//______________________________________________________________________________
TRotMatrix::TRotMatrix()
{
   // RotMatrix default constructor.
}


//______________________________________________________________________________
TRotMatrix::TRotMatrix(const char *name, const char *title, Double_t *matrix)
           :TNamed(name,title)
{
   // RotMatrix normal constructor.

   if (!matrix) { Error("ctor","No rotation is supplied"); return; }

   SetMatrix(matrix);
   if (!gGeometry) gGeometry = new TGeometry();
   fNumber = gGeometry->GetListOfMatrices()->GetSize();
   gGeometry->GetListOfMatrices()->Add(this);
}


//______________________________________________________________________________
TRotMatrix::TRotMatrix(const char *name, const char *title, Double_t theta, Double_t phi, Double_t psi)
           :TNamed(name,title)
{
   // RotMatrix normal constructor.

   printf("ERROR: This form of TRotMatrix constructor not implemented yet\n");

   Int_t i;
   fTheta  = theta;
   fPhi    = phi;
   fPsi    = psi;
   fType   = 2;
   for (i=0;i<9;i++) fMatrix[i] = 0;
   fMatrix[0] = 1;   fMatrix[4] = 1;   fMatrix[8] = 1;

   if (!gGeometry) gGeometry = new TGeometry();
   fNumber = gGeometry->GetListOfMatrices()->GetSize();
   gGeometry->GetListOfMatrices()->Add(this);
}


//______________________________________________________________________________
TRotMatrix::TRotMatrix(const char *name, const char *title, Double_t theta1, Double_t phi1
                                                  , Double_t theta2, Double_t phi2
                                                  , Double_t theta3, Double_t phi3)
                :TNamed(name,title)
{
   // RotMatrix normal constructor defined a la GEANT.
   //
   // The TRotMatrix constructor with six angles uses the GEANT convention:
   //
   // theta1 is the polar angle of the x-prim axis in the main reference system
   // (MRS), theta2 and theta3 have the same meaning for the y-prim and z-prim
   // axis.
   //
   // Phi1 is the azimuthal angle of the x-prim in the MRS and phi2 and phi3
   // have the same meaning for y-prim and z-prim.
   //
   //
   // for example, the unit matrix is defined in the following way.
   //
   //     x-prim || x, y-prim || y, z-prim || z
   //
   //     means:  theta1=90, theta2=90, theta3=0, phi1=0, phi2=90, phi3=0

   SetAngles(theta1,phi1,theta2,phi2,theta3,phi3);

   if (!gGeometry) gGeometry = new TGeometry();
   fNumber = gGeometry->GetListOfMatrices()->GetSize();
   gGeometry->GetListOfMatrices()->Add(this);
}


//______________________________________________________________________________
TRotMatrix::~TRotMatrix()
{
   // RotMatrix default destructor.

   if (gGeometry) gGeometry->GetListOfMatrices()->Remove(this);
}


//______________________________________________________________________________
Double_t  TRotMatrix::Determinant() const
{
   // Determinant() returns the value of the determiant of this matrix

   return
      fMatrix[0] * (fMatrix[4]*fMatrix[8] - fMatrix[7]*fMatrix[5])
    - fMatrix[3] * (fMatrix[1]*fMatrix[8] - fMatrix[7]*fMatrix[2])
    + fMatrix[6] * (fMatrix[1]*fMatrix[5] - fMatrix[4]*fMatrix[2]);
}


//______________________________________________________________________________
Double_t* TRotMatrix::GetGLMatrix(Double_t *rGLMatrix) const
{
   //  Convert this matrix to the OpenGL [4x4]
   //
   //  [  fMatrix[0]   fMatrix[1]   fMatrix[2]    0  ]
   //  [  fMatrix[3]   fMatrix[4]   fMatrix[5]    0  ]
   //  [  fMatrix[6]   fMatrix[7]   fMatrix[8]    0  ]
   //  [     0             0           0          1  ]
   //
   //  Input:
   //  -----
   //  Double_t *rGLMatrix - pointer to Double_t 4x4 buffer array
   //
   //  Return:
   //  ------
   //  Double_t pointer to the input buffer

   Double_t *glmatrix = rGLMatrix;
   const Double_t *matrix   = fMatrix;
   if (rGLMatrix)
   {
      for (Int_t i=0;i<3;i++) {
         for (Int_t j=0;j<3;j++) memcpy(glmatrix,matrix,3*sizeof(Double_t));
         matrix   += 3;
         glmatrix += 3;
         *glmatrix = 0.0;
         glmatrix++;
      }
      for (Int_t j=0;j<3;j++) {
         *glmatrix = 0.0;
         glmatrix++;
      }
      *glmatrix = 1.0;
   }
   return rGLMatrix;
}


//______________________________________________________________________________
const Double_t* TRotMatrix::SetAngles(Double_t theta1, Double_t phi1,
                Double_t theta2, Double_t phi2,Double_t theta3, Double_t phi3)
{
   // theta1 is the polar angle of the x-prim axis in the main reference system
   // (MRS), theta2 and theta3 have the same meaning for the y-prim and z-prim
   // axis.
   //
   // Phi1 is the azimuthal angle of the x-prim in the MRS and phi2 and phi3
   // have the same meaning for y-prim and z-prim.
   //
   //
   // for example, the unit matrix is defined in the following way.
   //
   //     x-prim || x, y-prim || y, z-prim || z
   //
   //     means:  theta1=90, theta2=90, theta3=0, phi1=0, phi2=90, phi3=0

   const Double_t degrad = 0.0174532925199432958;

   fTheta  = theta1;
   fPhi    = phi1;
   fPsi    = theta2;

   fType   = 2;
   if (!strcmp(GetName(),"Identity")) fType = 0;

   fMatrix[0] = TMath::Sin(theta1*degrad)*TMath::Cos(phi1*degrad);
   fMatrix[1] = TMath::Sin(theta1*degrad)*TMath::Sin(phi1*degrad);
   fMatrix[2] = TMath::Cos(theta1*degrad);
   fMatrix[3] = TMath::Sin(theta2*degrad)*TMath::Cos(phi2*degrad);
   fMatrix[4] = TMath::Sin(theta2*degrad)*TMath::Sin(phi2*degrad);
   fMatrix[5] = TMath::Cos(theta2*degrad);
   fMatrix[6] = TMath::Sin(theta3*degrad)*TMath::Cos(phi3*degrad);
   fMatrix[7] = TMath::Sin(theta3*degrad)*TMath::Sin(phi3*degrad);
   fMatrix[8] = TMath::Cos(theta3*degrad);

   SetReflection();
   return fMatrix;
}


//______________________________________________________________________________
void TRotMatrix::SetMatrix(const Double_t *matrix)
{
   // copy predefined 3x3 matrix into TRotMatrix object

   fTheta  = 0;
   fPhi    = 0;
   fPsi    = 0;
   fType   = 0;
   if (!matrix) return;
   fType   = 2;
   memcpy(fMatrix,matrix,9*sizeof(Double_t));
   SetReflection();
}


//______________________________________________________________________________
void TRotMatrix::SetReflection()
{
   // SetReflection() -  checks whether the determinant of this
   //                    matrix defines the reflection transformation
   //                    and set the "reflection" flag if any

   ResetBit(kReflection);
   if (Determinant() < 0) { fType=1; SetBit(kReflection);}
}


//______________________________________________________________________________
void TRotMatrix::Streamer(TBuffer &R__b)
{
   // Stream an object of class TRotMatrix.

   if (R__b.IsReading()) {
      UInt_t R__s, R__c;
      Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
      if (R__v > 1) {
         R__b.ReadClassBuffer(TRotMatrix::Class(), this, R__v, R__s, R__c);
         return;
      }
      //====process old versions before automatic schema evolution
      TNamed::Streamer(R__b);
      R__b >> fNumber;
      R__b >> fType;
      R__b >> fTheta;
      R__b >> fPhi;
      R__b >> fPsi;
      R__b.ReadStaticArray(fMatrix);
      R__b.CheckByteCount(R__s, R__c, TRotMatrix::IsA());
      //====end of old versions
      
   } else {
      R__b.WriteClassBuffer(TRotMatrix::Class(),this);
   }
}