//Begin_Html
/*
<img src="gif/t_transf.jpg">
*/
//End_Html
#include "Riostream.h"
#include "TObjArray.h"
#include "TGeoManager.h"
#include "TGeoMatrix.h"
TGeoIdentity *gGeoIdentity = 0;
const Int_t kN3 = 3*sizeof(Double_t);
const Int_t kN9 = 9*sizeof(Double_t);
ClassImp(TGeoMatrix)
TGeoMatrix::TGeoMatrix()
{
}
TGeoMatrix::TGeoMatrix(const TGeoMatrix &other)
:TNamed(other)
{
ResetBit(kGeoRegistered);
}
TGeoMatrix::TGeoMatrix(const char *name)
:TNamed(name, "")
{
}
TGeoMatrix::~TGeoMatrix()
{
if (IsRegistered() && gGeoManager) {
if (gGeoManager->GetListOfVolumes()) {
gGeoManager->GetListOfMatrices()->Remove(this);
Warning("dtor", "Registered matrix %s was removed", GetName());
}
}
}
TGeoMatrix& TGeoMatrix::operator = (const TGeoMatrix &matrix)
{
if (&matrix == this) return *this;
TNamed::operator=(matrix);
ResetBit(kGeoRegistered);
return *this;
}
TGeoMatrix &TGeoMatrix::operator*(const TGeoMatrix &right) const
{
static TGeoHMatrix h;
h = *this;
h.Multiply(&right);
return h;
}
Bool_t TGeoMatrix::operator ==(const TGeoMatrix &other) const
{
if (&other == this) return kTRUE;
Int_t i;
Bool_t tr1 = IsTranslation();
Bool_t tr2 = other.IsTranslation();
if ((tr1 & !tr2) || (tr2 & !tr1)) return kFALSE;
Bool_t rr1 = IsRotation();
Bool_t rr2 = other.IsRotation();
if ((rr1 & !rr2) || (rr2 & !rr1)) return kFALSE;
if (tr1) {
const Double_t *tr = GetTranslation();
const Double_t *otr = other.GetTranslation();
for (i=0; i<3; i++) if (TMath::Abs(tr[i]-otr[i])>1.E-10) return kFALSE;
}
if (rr1) {
const Double_t *rot = GetRotationMatrix();
const Double_t *orot = other.GetRotationMatrix();
for (i=0; i<9; i++) if (TMath::Abs(rot[i]-orot[i])>1.E-10) return kFALSE;
}
return kTRUE;
}
Bool_t TGeoMatrix::IsRotAboutZ() const
{
if (IsIdentity()) return kTRUE;
const Double_t *rot = GetRotationMatrix();
if (TMath::Abs(rot[6])>1E-9) return kFALSE;
if (TMath::Abs(rot[7])>1E-9) return kFALSE;
if ((1.-TMath::Abs(rot[8]))>1E-9) return kFALSE;
return kTRUE;
}
Int_t TGeoMatrix::GetByteCount() const
{
Int_t count = 4+28+strlen(GetName())+strlen(GetTitle());
if (IsTranslation()) count += 12;
if (IsScale()) count += 12;
if (IsCombi() || IsGeneral()) count += 4 + 36;
return count;
}
char *TGeoMatrix::GetPointerName() const
{
static char name[20];
sprintf(name,"pMatrix%d", GetUniqueID());
return name;
}
void TGeoMatrix::GetHomogenousMatrix(Double_t *hmat) const
{
Double_t *hmatrix = hmat;
const Double_t *mat = GetRotationMatrix();
for (Int_t i=0; i<3; i++) {
memcpy(hmatrix, mat, kN3);
mat += 3;
hmatrix += 3;
*hmatrix = 0.0;
hmatrix++;
}
memcpy(hmatrix, GetTranslation(), kN3);
hmatrix = hmat;
if (IsScale()) {
for (Int_t i=0; i<3; i++) {
*hmatrix *= GetScale()[i];
hmatrix += 5;
}
}
}
void TGeoMatrix::LocalToMaster(const Double_t *local, Double_t *master) const
{
if (IsIdentity()) {
memcpy(master, local, kN3);
return;
}
Int_t i;
const Double_t *tr = GetTranslation();
if (!IsRotation()) {
for (i=0; i<3; i++) master[i] = tr[i] + local[i];
return;
}
const Double_t *rot = GetRotationMatrix();
for (i=0; i<3; i++) {
master[i] = tr[i]
+ local[0]*rot[3*i]
+ local[1]*rot[3*i+1]
+ local[2]*rot[3*i+2];
}
}
void TGeoMatrix::LocalToMasterVect(const Double_t *local, Double_t *master) const
{
if (!IsRotation()) {
memcpy(master, local, kN3);
return;
}
const Double_t *rot = GetRotationMatrix();
for (Int_t i=0; i<3; i++) {
master[i] = local[0]*rot[3*i]
+ local[1]*rot[3*i+1]
+ local[2]*rot[3*i+2];
}
}
void TGeoMatrix::LocalToMasterBomb(const Double_t *local, Double_t *master) const
{
if (IsIdentity()) {
memcpy(master, local, kN3);
return;
}
Int_t i;
const Double_t *tr = GetTranslation();
Double_t bombtr[3];
gGeoManager->BombTranslation(tr, &bombtr[0]);
if (!IsRotation()) {
for (i=0; i<3; i++) master[i] = bombtr[i] + local[i];
return;
}
const Double_t *rot = GetRotationMatrix();
for (i=0; i<3; i++) {
master[i] = bombtr[i]
+ local[0]*rot[3*i]
+ local[1]*rot[3*i+1]
+ local[2]*rot[3*i+2];
}
}
void TGeoMatrix::MasterToLocal(const Double_t *master, Double_t *local) const
{
if (IsIdentity()) {
memcpy(local, master, kN3);
return;
}
const Double_t *tr = GetTranslation();
Double_t mt0 = master[0]-tr[0];
Double_t mt1 = master[1]-tr[1];
Double_t mt2 = master[2]-tr[2];
if (!IsRotation()) {
local[0] = mt0;
local[1] = mt1;
local[2] = mt2;
return;
}
const Double_t *rot = GetRotationMatrix();
local[0] = mt0*rot[0] + mt1*rot[3] + mt2*rot[6];
local[1] = mt0*rot[1] + mt1*rot[4] + mt2*rot[7];
local[2] = mt0*rot[2] + mt1*rot[5] + mt2*rot[8];
}
void TGeoMatrix::MasterToLocalVect(const Double_t *master, Double_t *local) const
{
if (!IsRotation()) {
memcpy(local, master, kN3);
return;
}
const Double_t *rot = GetRotationMatrix();
for (Int_t i=0; i<3; i++) {
local[i] = master[0]*rot[i]
+ master[1]*rot[i+3]
+ master[2]*rot[i+6];
}
}
void TGeoMatrix::MasterToLocalBomb(const Double_t *master, Double_t *local) const
{
if (IsIdentity()) {
memcpy(local, master, kN3);
return;
}
const Double_t *tr = GetTranslation();
Double_t bombtr[3];
Int_t i;
gGeoManager->UnbombTranslation(tr, &bombtr[0]);
if (!IsRotation()) {
for (i=0; i<3; i++) local[i] = master[i]-bombtr[i];
return;
}
const Double_t *rot = GetRotationMatrix();
for (i=0; i<3; i++) {
local[i] = (master[0]-bombtr[0])*rot[i]
+ (master[1]-bombtr[1])*rot[i+3]
+ (master[2]-bombtr[2])*rot[i+6];
}
}
void TGeoMatrix::Normalize(Double_t *vect)
{
Double_t normfactor = vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2];
if (normfactor <= 1E-10) return;
normfactor = 1./TMath::Sqrt(normfactor);
vect[0] *= normfactor;
vect[1] *= normfactor;
vect[2] *= normfactor;
}
void TGeoMatrix::Print(Option_t *) const
{
const Double_t *rot = GetRotationMatrix();
const Double_t *tr = GetTranslation();
printf("matrix %s - tr=%d rot=%d refl=%d scl=%d\n", GetName(),(Int_t)IsTranslation(),
(Int_t)IsRotation(), (Int_t)IsReflection(), (Int_t)IsScale());
printf("%10.6f%12.6f%12.6f Tx = %10.6f\n", rot[0], rot[1], rot[2], tr[0]);
printf("%10.6f%12.6f%12.6f Ty = %10.6f\n", rot[3], rot[4], rot[5], tr[1]);
printf("%10.6f%12.6f%12.6f Tz = %10.6f\n", rot[6], rot[7], rot[8], tr[2]);
if (IsScale()) {
const Double_t *scl = GetScale();
printf("Sx=%10.6fSy=%12.6fSz=%12.6f\n", scl[0], scl[1], scl[2]);
}
}
void TGeoMatrix::ReflectX(Bool_t, Bool_t)
{
}
void TGeoMatrix::ReflectY(Bool_t, Bool_t)
{
}
void TGeoMatrix::ReflectZ(Bool_t, Bool_t)
{
}
void TGeoMatrix::RegisterYourself()
{
if (!gGeoManager) {
Warning("RegisterYourself", "cannot register without geometry");
return;
}
if (!IsRegistered()) {
gGeoManager->RegisterMatrix(this);
SetBit(kGeoRegistered);
}
}
void TGeoMatrix::SetDefaultName()
{
if (!gGeoManager) return;
if (strlen(GetName())) return;
char type = 'n';
if (IsTranslation()) type = 't';
if (IsRotation()) type = 'r';
if (IsScale()) type = 's';
if (IsCombi()) type = 'c';
if (IsGeneral()) type = 'g';
TObjArray *matrices = gGeoManager->GetListOfMatrices();
Int_t index = 0;
if (matrices) index =matrices->GetEntriesFast() - 1;
Int_t digits = 1;
Int_t num = 10;
while ((Int_t)(index/num)) {
digits++;
num *= 10;
}
char *name = new char[digits+2];
sprintf(name, "%c%i", type, index);
SetName(name);
}
ClassImp(TGeoTranslation)
TGeoTranslation::TGeoTranslation()
{
for (Int_t i=0; i<3; i++) fTranslation[i] = 0;
}
TGeoTranslation::TGeoTranslation(const TGeoTranslation &other)
:TGeoMatrix(other)
{
SetTranslation(other);
}
TGeoTranslation::TGeoTranslation(const TGeoMatrix &other)
:TGeoMatrix(other)
{
SetTranslation(other);
}
TGeoTranslation::TGeoTranslation(Double_t dx, Double_t dy, Double_t dz)
:TGeoMatrix("")
{
if (dx || dy || dz) SetBit(kGeoTranslation);
SetTranslation(dx, dy, dz);
}
TGeoTranslation::TGeoTranslation(const char *name, Double_t dx, Double_t dy, Double_t dz)
:TGeoMatrix(name)
{
if (dx || dy || dz) SetBit(kGeoTranslation);
SetTranslation(dx, dy, dz);
}
TGeoTranslation& TGeoTranslation::operator = (const TGeoMatrix &matrix)
{
if (&matrix == this) return *this;
TGeoMatrix::operator=(matrix);
SetTranslation(matrix);
return *this;
}
TGeoMatrix& TGeoTranslation::Inverse() const
{
static TGeoHMatrix h;
h = *this;
Double_t tr[3];
tr[0] = -fTranslation[0];
tr[1] = -fTranslation[1];
tr[2] = -fTranslation[2];
h.SetTranslation(tr);
return h;
}
void TGeoTranslation::Add(const TGeoTranslation *other)
{
const Double_t *trans = other->GetTranslation();
for (Int_t i=0; i<3; i++)
fTranslation[i] += trans[i];
}
TGeoMatrix *TGeoTranslation::MakeClone() const
{
TGeoMatrix *matrix = new TGeoTranslation(*this);
return matrix;
}
void TGeoTranslation::RotateX(Double_t )
{
Warning("RotateX", "Not implemented. Use TGeoCombiTrans instead");
}
void TGeoTranslation::RotateY(Double_t )
{
Warning("RotateY", "Not implemented. Use TGeoCombiTrans instead");
}
void TGeoTranslation::RotateZ(Double_t )
{
Warning("RotateZ", "Not implemented. Use TGeoCombiTrans instead");
}
void TGeoTranslation::Subtract(const TGeoTranslation *other)
{
const Double_t *trans = other->GetTranslation();
for (Int_t i=0; i<3; i++)
fTranslation[i] -= trans[i];
}
void TGeoTranslation::SetTranslation(Double_t dx, Double_t dy, Double_t dz)
{
fTranslation[0] = dx;
fTranslation[1] = dy;
fTranslation[2] = dz;
if (dx || dy || dz) SetBit(kGeoTranslation);
else ResetBit(kGeoTranslation);
}
void TGeoTranslation::SetTranslation(const TGeoMatrix &other)
{
SetBit(kGeoTranslation, other.IsTranslation());
const Double_t *transl = other.GetTranslation();
memcpy(fTranslation, transl, kN3);
}
void TGeoTranslation::LocalToMaster(const Double_t *local, Double_t *master) const
{
const Double_t *tr = GetTranslation();
for (Int_t i=0; i<3; i++)
master[i] = tr[i] + local[i];
}
void TGeoTranslation::LocalToMasterVect(const Double_t *local, Double_t *master) const
{
memcpy(master, local, kN3);
}
void TGeoTranslation::LocalToMasterBomb(const Double_t *local, Double_t *master) const
{
const Double_t *tr = GetTranslation();
Double_t bombtr[3];
gGeoManager->BombTranslation(tr, &bombtr[0]);
for (Int_t i=0; i<3; i++)
master[i] = bombtr[i] + local[i];
}
void TGeoTranslation::MasterToLocal(const Double_t *master, Double_t *local) const
{
const Double_t *tr = GetTranslation();
for (Int_t i=0; i<3; i++)
local[i] = master[i]-tr[i];
}
void TGeoTranslation::MasterToLocalVect(const Double_t *master, Double_t *local) const
{
memcpy(local, master, kN3);
}
void TGeoTranslation::MasterToLocalBomb(const Double_t *master, Double_t *local) const
{
const Double_t *tr = GetTranslation();
Double_t bombtr[3];
gGeoManager->UnbombTranslation(tr, &bombtr[0]);
for (Int_t i=0; i<3; i++)
local[i] = master[i]-bombtr[i];
}
void TGeoTranslation::SavePrimitive(ostream &out, Option_t * )
{
if (TestBit(kGeoSavePrimitive)) return;
out << " // Translation: " << GetName() << endl;
out << " dx = " << fTranslation[0] << ";" << endl;
out << " dy = " << fTranslation[1] << ";" << endl;
out << " dz = " << fTranslation[2] << ";" << endl;
out << " TGeoTranslation *" << GetPointerName() << " = new TGeoTranslation(\"" << GetName() << "\",dx,dy,dz);" << endl;
TObject::SetBit(kGeoSavePrimitive);
}
ClassImp(TGeoRotation)
TGeoRotation::TGeoRotation()
{
for (Int_t i=0; i<9; i++) {
if (i%4) fRotationMatrix[i] = 0;
else fRotationMatrix[i] = 1.0;
}
}
TGeoRotation::TGeoRotation(const TGeoRotation &other)
:TGeoMatrix(other)
{
SetRotation(other);
}
TGeoRotation::TGeoRotation(const TGeoMatrix &other)
:TGeoMatrix(other)
{
SetRotation(other);
}
TGeoRotation::TGeoRotation(const char *name)
:TGeoMatrix(name)
{
for (Int_t i=0; i<9; i++) {
if (i%4) fRotationMatrix[i] = 0;
else fRotationMatrix[i] = 1.0;
}
}
TGeoRotation::TGeoRotation(const char *name, Double_t phi, Double_t theta, Double_t psi)
:TGeoMatrix(name)
{
SetAngles(phi, theta, psi);
}
TGeoRotation::TGeoRotation(const char *name, Double_t theta1, Double_t phi1, Double_t theta2, Double_t phi2,
Double_t theta3, Double_t phi3)
:TGeoMatrix(name)
{
SetAngles(theta1, phi1, theta2, phi2, theta3, phi3);
}
TGeoRotation& TGeoRotation::operator = (const TGeoMatrix &other)
{
if (&other == this) return *this;
TGeoMatrix::operator=(other);
SetRotation(other);
return *this;
}
TGeoMatrix& TGeoRotation::Inverse() const
{
static TGeoHMatrix h;
h = *this;
Double_t newrot[9];
newrot[0] = fRotationMatrix[0];
newrot[1] = fRotationMatrix[3];
newrot[2] = fRotationMatrix[6];
newrot[3] = fRotationMatrix[1];
newrot[4] = fRotationMatrix[4];
newrot[5] = fRotationMatrix[7];
newrot[6] = fRotationMatrix[2];
newrot[7] = fRotationMatrix[5];
newrot[8] = fRotationMatrix[8];
h.SetRotation(newrot);
return h;
}
Bool_t TGeoRotation::IsValid() const
{
const Double_t *r = fRotationMatrix;
Double_t cij;
for (Int_t i=0; i<2; i++) {
for (Int_t j=i+1; j<3; j++) {
cij = TMath::Abs(r[i]*r[j]+r[i+3]*r[j+3]+r[i+6]*r[j+6]);
if (cij>1E-4) return kFALSE;
cij = TMath::Abs(r[3*i]*r[3*j]+r[3*i+1]*r[3*j+1]+r[3*i+2]*r[3*j+2]);
if (cij>1E-4) return kFALSE;
}
}
return kTRUE;
}
void TGeoRotation::Clear(Option_t *)
{
memcpy(fRotationMatrix,kIdentityMatrix,kN9);
ResetBit(kGeoRotation);
}
void TGeoRotation::FastRotZ(Double_t *sincos)
{
fRotationMatrix[0] = sincos[1];
fRotationMatrix[1] = -sincos[0];
fRotationMatrix[3] = sincos[0];
fRotationMatrix[4] = sincos[1];
SetBit(kGeoRotation);
}
Double_t TGeoRotation::GetPhiRotation(Bool_t fixX) const
{
Double_t phi;
if (fixX) phi = 180.*TMath::ATan2(-fRotationMatrix[1],fRotationMatrix[4])/TMath::Pi();
else phi = 180.*TMath::ATan2(fRotationMatrix[3], fRotationMatrix[0])/TMath::Pi();
return phi;
}
void TGeoRotation::LocalToMaster(const Double_t *local, Double_t *master) const
{
const Double_t *rot = GetRotationMatrix();
for (Int_t i=0; i<3; i++) {
master[i] = local[0]*rot[3*i]
+ local[1]*rot[3*i+1]
+ local[2]*rot[3*i+2];
}
}
void TGeoRotation::MasterToLocal(const Double_t *master, Double_t *local) const
{
const Double_t *rot = GetRotationMatrix();
for (Int_t i=0; i<3; i++) {
local[i] = master[0]*rot[i]
+ master[1]*rot[i+3]
+ master[2]*rot[i+6];
}
}
TGeoMatrix *TGeoRotation::MakeClone() const
{
TGeoMatrix *matrix = new TGeoRotation(*this);
return matrix;
}
void TGeoRotation::RotateX(Double_t angle)
{
SetBit(kGeoRotation);
Double_t phi = angle*TMath::DegToRad();
Double_t c = TMath::Cos(phi);
Double_t s = TMath::Sin(phi);
Double_t v[9];
v[0] = fRotationMatrix[0];
v[1] = fRotationMatrix[1];
v[2] = fRotationMatrix[2];
v[3] = c*fRotationMatrix[3]-s*fRotationMatrix[6];
v[4] = c*fRotationMatrix[4]-s*fRotationMatrix[7];
v[5] = c*fRotationMatrix[5]-s*fRotationMatrix[8];
v[6] = s*fRotationMatrix[3]+c*fRotationMatrix[6];
v[7] = s*fRotationMatrix[4]+c*fRotationMatrix[7];
v[8] = s*fRotationMatrix[5]+c*fRotationMatrix[8];
memcpy(fRotationMatrix, v, kN9);
}
void TGeoRotation::RotateY(Double_t angle)
{
SetBit(kGeoRotation);
Double_t phi = angle*TMath::DegToRad();
Double_t c = TMath::Cos(phi);
Double_t s = TMath::Sin(phi);
Double_t v[9];
v[0] = c*fRotationMatrix[0]+s*fRotationMatrix[6];
v[1] = c*fRotationMatrix[1]+s*fRotationMatrix[7];
v[2] = c*fRotationMatrix[2]+s*fRotationMatrix[8];
v[3] = fRotationMatrix[3];
v[4] = fRotationMatrix[4];
v[5] = fRotationMatrix[5];
v[6] = -s*fRotationMatrix[0]+c*fRotationMatrix[6];
v[7] = -s*fRotationMatrix[1]+c*fRotationMatrix[7];
v[8] = -s*fRotationMatrix[2]+c*fRotationMatrix[8];
memcpy(fRotationMatrix, v, kN9);
}
void TGeoRotation::RotateZ(Double_t angle)
{
SetBit(kGeoRotation);
Double_t phi = angle*TMath::DegToRad();
Double_t c = TMath::Cos(phi);
Double_t s = TMath::Sin(phi);
Double_t v[9];
v[0] = c*fRotationMatrix[0]-s*fRotationMatrix[3];
v[1] = c*fRotationMatrix[1]-s*fRotationMatrix[4];
v[2] = c*fRotationMatrix[2]-s*fRotationMatrix[5];
v[3] = s*fRotationMatrix[0]+c*fRotationMatrix[3];
v[4] = s*fRotationMatrix[1]+c*fRotationMatrix[4];
v[5] = s*fRotationMatrix[2]+c*fRotationMatrix[5];
v[6] = fRotationMatrix[6];
v[7] = fRotationMatrix[7];
v[8] = fRotationMatrix[8];
memcpy(&fRotationMatrix[0],v,kN9);
}
void TGeoRotation::ReflectX(Bool_t leftside, Bool_t)
{
if (leftside) {
fRotationMatrix[0]=-fRotationMatrix[0];
fRotationMatrix[1]=-fRotationMatrix[1];
fRotationMatrix[2]=-fRotationMatrix[2];
} else {
fRotationMatrix[0]=-fRotationMatrix[0];
fRotationMatrix[3]=-fRotationMatrix[3];
fRotationMatrix[6]=-fRotationMatrix[6];
}
SetBit(kGeoRotation);
SetBit(kGeoReflection, !IsReflection());
}
void TGeoRotation::ReflectY(Bool_t leftside, Bool_t)
{
if (leftside) {
fRotationMatrix[3]=-fRotationMatrix[3];
fRotationMatrix[4]=-fRotationMatrix[4];
fRotationMatrix[5]=-fRotationMatrix[5];
} else {
fRotationMatrix[1]=-fRotationMatrix[1];
fRotationMatrix[4]=-fRotationMatrix[4];
fRotationMatrix[7]=-fRotationMatrix[7];
}
SetBit(kGeoRotation);
SetBit(kGeoReflection, !IsReflection());
}
void TGeoRotation::ReflectZ(Bool_t leftside, Bool_t)
{
if (leftside) {
fRotationMatrix[6]=-fRotationMatrix[6];
fRotationMatrix[7]=-fRotationMatrix[7];
fRotationMatrix[8]=-fRotationMatrix[8];
} else {
fRotationMatrix[2]=-fRotationMatrix[2];
fRotationMatrix[5]=-fRotationMatrix[5];
fRotationMatrix[8]=-fRotationMatrix[8];
}
SetBit(kGeoRotation);
SetBit(kGeoReflection, !IsReflection());
}
void TGeoRotation::SavePrimitive(ostream &out, Option_t * )
{
if (TestBit(kGeoSavePrimitive)) return;
out << " // Rotation: " << GetName() << endl;
Double_t th1,ph1,th2,ph2,th3,ph3;
GetAngles(th1,ph1,th2,ph2,th3,ph3);
out << " thx = " << th1 << "; phx = " << ph1 << ";" << endl;
out << " thy = " << th2 << "; phy = " << ph2 << ";" << endl;
out << " thz = " << th3 << "; phz = " << ph3 << ";" << endl;
out << " TGeoRotation *" << GetPointerName() << " = new TGeoRotation(\"" << GetName() << "\",thx,phx,thy,phy,thz,phz);" << endl;
TObject::SetBit(kGeoSavePrimitive);
}
void TGeoRotation::SetRotation(const TGeoMatrix &other)
{
SetBit(kGeoRotation, other.IsRotation());
const Double_t *rot = other.GetRotationMatrix();
SetMatrix(rot);
}
void TGeoRotation::SetAngles(Double_t phi, Double_t theta, Double_t psi)
{
Double_t degrad = TMath::Pi()/180.;
Double_t sinphi = TMath::Sin(degrad*phi);
Double_t cosphi = TMath::Cos(degrad*phi);
Double_t sinthe = TMath::Sin(degrad*theta);
Double_t costhe = TMath::Cos(degrad*theta);
Double_t sinpsi = TMath::Sin(degrad*psi);
Double_t cospsi = TMath::Cos(degrad*psi);
fRotationMatrix[0] = cospsi*cosphi - costhe*sinphi*sinpsi;
fRotationMatrix[1] = -sinpsi*cosphi - costhe*sinphi*cospsi;
fRotationMatrix[2] = sinthe*sinphi;
fRotationMatrix[3] = cospsi*sinphi + costhe*cosphi*sinpsi;
fRotationMatrix[4] = -sinpsi*sinphi + costhe*cosphi*cospsi;
fRotationMatrix[5] = -sinthe*cosphi;
fRotationMatrix[6] = sinpsi*sinthe;
fRotationMatrix[7] = cospsi*sinthe;
fRotationMatrix[8] = costhe;
if (!IsValid()) Error("SetAngles", "invalid rotation (Euler angles : phi=%f theta=%f psi=%f)",phi,theta,psi);
CheckMatrix();
}
void TGeoRotation::SetAngles(Double_t theta1, Double_t phi1, Double_t theta2, Double_t phi2,
Double_t theta3, Double_t phi3)
{
Double_t degrad = TMath::Pi()/180.;
fRotationMatrix[0] = TMath::Cos(degrad*phi1)*TMath::Sin(degrad*theta1);
fRotationMatrix[3] = TMath::Sin(degrad*phi1)*TMath::Sin(degrad*theta1);
fRotationMatrix[6] = TMath::Cos(degrad*theta1);
fRotationMatrix[1] = TMath::Cos(degrad*phi2)*TMath::Sin(degrad*theta2);
fRotationMatrix[4] = TMath::Sin(degrad*phi2)*TMath::Sin(degrad*theta2);
fRotationMatrix[7] = TMath::Cos(degrad*theta2);
fRotationMatrix[2] = TMath::Cos(degrad*phi3)*TMath::Sin(degrad*theta3);
fRotationMatrix[5] = TMath::Sin(degrad*phi3)*TMath::Sin(degrad*theta3);
fRotationMatrix[8] = TMath::Cos(degrad*theta3);
for (Int_t i=0; i<9; i++) {
if (TMath::Abs(fRotationMatrix[i])<1E-15) fRotationMatrix[i] = 0;
if (TMath::Abs(fRotationMatrix[i]-1)<1E-15) fRotationMatrix[i] = 1;
if (TMath::Abs(fRotationMatrix[i]+1)<1E-15) fRotationMatrix[i] = -1;
}
if (!IsValid()) Error("SetAngles", "invalid rotation (G3 angles, th1=%f phi1=%f, th2=%f ph2=%f, th3=%f phi3=%f)",
theta1,phi1,theta2,phi2,theta3,phi3);
CheckMatrix();
}
void TGeoRotation::GetAngles(Double_t &theta1, Double_t &phi1, Double_t &theta2, Double_t &phi2,
Double_t &theta3, Double_t &phi3) const
{
Double_t raddeg = 180./TMath::Pi();
theta1 = raddeg*TMath::ACos(fRotationMatrix[6]);
theta2 = raddeg*TMath::ACos(fRotationMatrix[7]);
theta3 = raddeg*TMath::ACos(fRotationMatrix[8]);
if (TMath::Abs(fRotationMatrix[0])<1E-6 && TMath::Abs(fRotationMatrix[3])<1E-6) phi1=0.;
else phi1 = raddeg*TMath::ATan2(fRotationMatrix[3],fRotationMatrix[0]);
if (phi1<0) phi1+=360.;
if (TMath::Abs(fRotationMatrix[1])<1E-6 && TMath::Abs(fRotationMatrix[4])<1E-6) phi2=0.;
else phi2 = raddeg*TMath::ATan2(fRotationMatrix[4],fRotationMatrix[1]);
if (phi2<0) phi2+=360.;
if (TMath::Abs(fRotationMatrix[2])<1E-6 && TMath::Abs(fRotationMatrix[5])<1E-6) phi3=0.;
else phi3 = raddeg*TMath::ATan2(fRotationMatrix[5],fRotationMatrix[2]);
if (phi3<0) phi3+=360.;
}
void TGeoRotation::GetAngles(Double_t &phi, Double_t &theta, Double_t &psi) const
{
const Double_t *m = fRotationMatrix;
if (TMath::Abs(1.-TMath::Abs(m[8]))<1.e-9) {
theta = TMath::ACos(m[8])*TMath::RadToDeg();
phi = TMath::ATan2(-m[8]*m[1],m[0])*TMath::RadToDeg();
psi = 0.;
return;
}
phi = TMath::ATan2(m[2],-m[5]);
Double_t sphi = TMath::Sin(phi);
if (TMath::Abs(sphi)<1.e-9) theta = -TMath::ASin(m[5]/TMath::Cos(phi))*TMath::RadToDeg();
else theta = TMath::ASin(m[2]/sphi)*TMath::RadToDeg();
phi *= TMath::RadToDeg();
psi = TMath::ATan2(m[6],m[7])*TMath::RadToDeg();
}
Double_t TGeoRotation::Determinant() const
{
Double_t
det = fRotationMatrix[0]*fRotationMatrix[4]*fRotationMatrix[8] +
fRotationMatrix[3]*fRotationMatrix[7]*fRotationMatrix[2] +
fRotationMatrix[6]*fRotationMatrix[1]*fRotationMatrix[5] -
fRotationMatrix[2]*fRotationMatrix[4]*fRotationMatrix[6] -
fRotationMatrix[5]*fRotationMatrix[7]*fRotationMatrix[0] -
fRotationMatrix[8]*fRotationMatrix[1]*fRotationMatrix[3];
return det;
}
void TGeoRotation::CheckMatrix()
{
if (Determinant() < 0) SetBit(kGeoReflection);
Double_t dd = fRotationMatrix[0] + fRotationMatrix[4] + fRotationMatrix[8] - 3.;
if (TMath::Abs(dd) < 1.E-12) ResetBit(kGeoRotation);
else SetBit(kGeoRotation);
}
void TGeoRotation::GetInverse(Double_t *invmat) const
{
if (!invmat) {
Error("GetInverse", "no place to store the inverse matrix");
}
for (Int_t i=0; i<3; i++) {
for (Int_t j=0; j<3; j++) {
invmat[3*i+j] = fRotationMatrix[3*j+i];
}
}
}
void TGeoRotation::MultiplyBy(TGeoRotation *rot, Bool_t after)
{
const Double_t *matleft, *matright;
SetBit(kGeoRotation);
Double_t newmat[9] = {0};
if (after) {
matleft = &fRotationMatrix[0];
matright = rot->GetRotationMatrix();
} else {
matleft = rot->GetRotationMatrix();
matright = &fRotationMatrix[0];
}
for (Int_t i=0; i<3; i++) {
for (Int_t j=0; j<3; j++) {
for (Int_t k=0; k<3; k++) {
newmat[3*i+j] += matleft[3*i+k] * matright[3*k+j];
}
}
}
memcpy(&fRotationMatrix[0], &newmat[0], kN9);
}
ClassImp(TGeoScale)
TGeoScale::TGeoScale()
{
SetBit(kGeoScale);
for (Int_t i=0; i<3; i++) fScale[i] = 1.;
}
TGeoScale::TGeoScale(const TGeoScale &other)
:TGeoMatrix(other)
{
SetBit(kGeoScale);
const Double_t *scl = other.GetScale();
memcpy(fScale, scl, kN3);
if (fScale[0]*fScale[1]*fScale[2]<0) SetBit(kGeoReflection);
else SetBit(kGeoReflection, kFALSE);
}
TGeoScale::TGeoScale(Double_t sx, Double_t sy, Double_t sz)
:TGeoMatrix("")
{
SetBit(kGeoScale);
SetScale(sx, sy, sz);
}
TGeoScale::TGeoScale(const char *name, Double_t sx, Double_t sy, Double_t sz)
:TGeoMatrix(name)
{
SetBit(kGeoScale);
SetScale(sx, sy, sz);
}
TGeoScale::~TGeoScale()
{
}
TGeoMatrix& TGeoScale::Inverse() const
{
static TGeoHMatrix h;
h = *this;
Double_t scale[3];
scale[0] = 1./fScale[0];
scale[1] = 1./fScale[1];
scale[2] = 1./fScale[2];
h.SetScale(scale);
return h;
}
void TGeoScale::SetScale(Double_t sx, Double_t sy, Double_t sz)
{
if (TMath::Abs(sx*sy*sz) < 1.E-10) {
Error("SetScale", "Invalid scale %f, %f, %f for transformation %s",sx,sy,sx,GetName());
return;
}
fScale[0] = sx;
fScale[1] = sy;
fScale[2] = sz;
if (sx*sy*sz<0) SetBit(kGeoReflection);
else SetBit(kGeoReflection, kFALSE);
}
void TGeoScale::LocalToMaster(const Double_t *local, Double_t *master) const
{
master[0] = local[0]*fScale[0];
master[1] = local[1]*fScale[1];
master[2] = local[2]*fScale[2];
}
Double_t TGeoScale::LocalToMaster(Double_t dist, const Double_t *dir) const
{
Double_t scale;
if (!dir) {
scale = TMath::Abs(fScale[0]);
if (TMath::Abs(fScale[1])<scale) scale = TMath::Abs(fScale[1]);
if (TMath::Abs(fScale[2])<scale) scale = TMath::Abs(fScale[2]);
} else {
scale = fScale[0]*fScale[0]*dir[0]*dir[0] +
fScale[1]*fScale[1]*dir[1]*dir[1] +
fScale[2]*fScale[2]*dir[2]*dir[2];
scale = TMath::Sqrt(scale);
}
return scale*dist;
}
TGeoMatrix *TGeoScale::MakeClone() const
{
TGeoMatrix *matrix = new TGeoScale(*this);
return matrix;
}
void TGeoScale::MasterToLocal(const Double_t *master, Double_t *local) const
{
local[0] = master[0]/fScale[0];
local[1] = master[1]/fScale[1];
local[2] = master[2]/fScale[2];
}
Double_t TGeoScale::MasterToLocal(Double_t dist, const Double_t *dir) const
{
Double_t scale;
if (!dir) {
scale = TMath::Abs(fScale[0]);
if (TMath::Abs(fScale[1])>scale) scale = TMath::Abs(fScale[1]);
if (TMath::Abs(fScale[2])>scale) scale = TMath::Abs(fScale[2]);
scale = 1./scale;
} else {
scale = (dir[0]*dir[0])/(fScale[0]*fScale[0]) +
(dir[1]*dir[1])/(fScale[1]*fScale[1]) +
(dir[2]*dir[2])/(fScale[2]*fScale[2]);
scale = TMath::Sqrt(scale);
}
return scale*dist;
}
ClassImp(TGeoCombiTrans)
TGeoCombiTrans::TGeoCombiTrans()
{
for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
fRotation = 0;
}
TGeoCombiTrans::TGeoCombiTrans(const TGeoCombiTrans &other)
:TGeoMatrix(other)
{
Int_t i;
if (other.IsTranslation()) {
const Double_t *trans = other.GetTranslation();
memcpy(fTranslation, trans, kN3);
} else {
for (i=0; i<3; i++) fTranslation[i] = 0.0;
}
if (other.IsRotation()) {
const TGeoRotation rot = *other.GetRotation();
fRotation = new TGeoRotation(rot);
SetBit(kGeoMatrixOwned);
}
else fRotation = 0;
}
TGeoCombiTrans::TGeoCombiTrans(const TGeoMatrix &other)
:TGeoMatrix(other)
{
Int_t i;
if (other.IsTranslation()) {
SetBit(kGeoTranslation);
memcpy(fTranslation,other.GetTranslation(),kN3);
} else {
for (i=0; i<3; i++) fTranslation[i] = 0.0;
}
if (other.IsRotation()) {
SetBit(kGeoRotation);
SetBit(kGeoMatrixOwned);
fRotation = new TGeoRotation(other);
}
else fRotation = 0;
}
TGeoCombiTrans::TGeoCombiTrans(const TGeoTranslation &tr, const TGeoRotation &rot)
{
if (tr.IsTranslation()) {
SetBit(kGeoTranslation);
const Double_t *trans = tr.GetTranslation();
memcpy(fTranslation, trans, kN3);
} else {
for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
}
if (rot.IsRotation()) {
SetBit(kGeoRotation);
SetBit(kGeoMatrixOwned);
fRotation = new TGeoRotation(rot);
SetBit(kGeoReflection, rot.TestBit(kGeoReflection));
}
else fRotation = 0;
}
TGeoCombiTrans::TGeoCombiTrans(const char *name)
:TGeoMatrix(name)
{
for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
fRotation = 0;
}
TGeoCombiTrans::TGeoCombiTrans(Double_t dx, Double_t dy, Double_t dz, TGeoRotation *rot)
:TGeoMatrix("")
{
SetTranslation(dx, dy, dz);
fRotation = 0;
SetRotation(rot);
}
TGeoCombiTrans::TGeoCombiTrans(const char *name, Double_t dx, Double_t dy, Double_t dz, TGeoRotation *rot)
:TGeoMatrix(name)
{
SetTranslation(dx, dy, dz);
fRotation = 0;
SetRotation(rot);
}
TGeoCombiTrans &TGeoCombiTrans::operator=(const TGeoMatrix &matrix)
{
if (&matrix == this) return *this;
Clear();
TGeoMatrix::operator=(matrix);
if (matrix.IsTranslation()) {
SetBit(kGeoTranslation);
memcpy(fTranslation,matrix.GetTranslation(),kN3);
}
if (matrix.IsRotation()) {
SetBit(kGeoRotation);
if (!fRotation) {
fRotation = new TGeoRotation();
SetBit(kGeoMatrixOwned);
} else {
if (!TestBit(kGeoMatrixOwned)) {
fRotation = new TGeoRotation();
SetBit(kGeoMatrixOwned);
}
}
fRotation->SetMatrix(matrix.GetRotationMatrix());
fRotation->SetBit(kGeoReflection, matrix.TestBit(kGeoReflection));
fRotation->SetBit(kGeoRotation);
} else {
if (fRotation && TestBit(kGeoMatrixOwned)) delete fRotation;
ResetBit(kGeoMatrixOwned);
fRotation = 0;
}
return *this;
}
TGeoCombiTrans::~TGeoCombiTrans()
{
if (fRotation) {
if(TestBit(TGeoMatrix::kGeoMatrixOwned) && !fRotation->IsRegistered()) delete fRotation;
}
}
void TGeoCombiTrans::Clear(Option_t *)
{
if (IsTranslation()) {
ResetBit(kGeoTranslation);
memset(fTranslation, 0, kN3);
}
if (fRotation) {
if (TestBit(kGeoMatrixOwned)) delete fRotation;
fRotation = 0;
}
ResetBit(kGeoRotation);
ResetBit(kGeoReflection);
ResetBit(kGeoMatrixOwned);
}
TGeoMatrix& TGeoCombiTrans::Inverse() const
{
static TGeoHMatrix h;
h = *this;
Bool_t is_tr = IsTranslation();
Bool_t is_rot = IsRotation();
Double_t tr[3];
Double_t newrot[9];
const Double_t *rot = GetRotationMatrix();
tr[0] = -fTranslation[0]*rot[0] - fTranslation[1]*rot[3] - fTranslation[2]*rot[6];
tr[1] = -fTranslation[0]*rot[1] - fTranslation[1]*rot[4] - fTranslation[2]*rot[7];
tr[2] = -fTranslation[0]*rot[2] - fTranslation[1]*rot[5] - fTranslation[2]*rot[8];
h.SetTranslation(tr);
newrot[0] = rot[0];
newrot[1] = rot[3];
newrot[2] = rot[6];
newrot[3] = rot[1];
newrot[4] = rot[4];
newrot[5] = rot[7];
newrot[6] = rot[2];
newrot[7] = rot[5];
newrot[8] = rot[8];
h.SetRotation(newrot);
h.SetBit(kGeoTranslation,is_tr);
h.SetBit(kGeoRotation,is_rot);
return h;
}
TGeoMatrix *TGeoCombiTrans::MakeClone() const
{
TGeoMatrix *matrix = new TGeoCombiTrans(*this);
return matrix;
}
void TGeoCombiTrans::RegisterYourself()
{
TGeoMatrix::RegisterYourself();
if (fRotation && fRotation->IsRotation()) fRotation->RegisterYourself();
}
void TGeoCombiTrans::RotateX(Double_t angle)
{
if (!fRotation || !TestBit(kGeoMatrixOwned)) {
if (fRotation) fRotation = new TGeoRotation(*fRotation);
else fRotation = new TGeoRotation();
SetBit(kGeoMatrixOwned);
}
SetBit(kGeoRotation);
const Double_t *rot = fRotation->GetRotationMatrix();
Double_t phi = angle*TMath::DegToRad();
Double_t c = TMath::Cos(phi);
Double_t s = TMath::Sin(phi);
Double_t v[9];
v[0] = rot[0];
v[1] = rot[1];
v[2] = rot[2];
v[3] = c*rot[3]-s*rot[6];
v[4] = c*rot[4]-s*rot[7];
v[5] = c*rot[5]-s*rot[8];
v[6] = s*rot[3]+c*rot[6];
v[7] = s*rot[4]+c*rot[7];
v[8] = s*rot[5]+c*rot[8];
fRotation->SetMatrix(v);
fRotation->SetBit(kGeoRotation);
if (!IsTranslation()) return;
v[0] = fTranslation[0];
v[1] = c*fTranslation[1]-s*fTranslation[2];
v[2] = s*fTranslation[1]+c*fTranslation[2];
memcpy(fTranslation,v,kN3);
}
void TGeoCombiTrans::RotateY(Double_t angle)
{
if (!fRotation || !TestBit(kGeoMatrixOwned)) {
if (fRotation) fRotation = new TGeoRotation(*fRotation);
else fRotation = new TGeoRotation();
SetBit(kGeoMatrixOwned);
}
SetBit(kGeoRotation);
const Double_t *rot = fRotation->GetRotationMatrix();
Double_t phi = angle*TMath::DegToRad();
Double_t c = TMath::Cos(phi);
Double_t s = TMath::Sin(phi);
Double_t v[9];
v[0] = c*rot[0]+s*rot[6];
v[1] = c*rot[1]+s*rot[7];
v[2] = c*rot[2]+s*rot[8];
v[3] = rot[3];
v[4] = rot[4];
v[5] = rot[5];
v[6] = -s*rot[0]+c*rot[6];
v[7] = -s*rot[1]+c*rot[7];
v[8] = -s*rot[2]+c*rot[8];
fRotation->SetMatrix(v);
fRotation->SetBit(kGeoRotation);
if (!IsTranslation()) return;
v[0] = c*fTranslation[0]+s*fTranslation[2];
v[1] = fTranslation[1];
v[2] = -s*fTranslation[0]+c*fTranslation[2];
memcpy(fTranslation,v,kN3);
}
void TGeoCombiTrans::RotateZ(Double_t angle)
{
if (!fRotation || !TestBit(kGeoMatrixOwned)) {
if (fRotation) fRotation = new TGeoRotation(*fRotation);
else fRotation = new TGeoRotation();
SetBit(kGeoMatrixOwned);
}
SetBit(kGeoRotation);
const Double_t *rot = fRotation->GetRotationMatrix();
Double_t phi = angle*TMath::DegToRad();
Double_t c = TMath::Cos(phi);
Double_t s = TMath::Sin(phi);
Double_t v[9];
v[0] = c*rot[0]-s*rot[3];
v[1] = c*rot[1]-s*rot[4];
v[2] = c*rot[2]-s*rot[5];
v[3] = s*rot[0]+c*rot[3];
v[4] = s*rot[1]+c*rot[4];
v[5] = s*rot[2]+c*rot[5];
v[6] = rot[6];
v[7] = rot[7];
v[8] = rot[8];
fRotation->SetMatrix(v);
fRotation->SetBit(kGeoRotation);
if (!IsTranslation()) return;
v[0] = c*fTranslation[0]-s*fTranslation[1];
v[1] = s*fTranslation[0]+c*fTranslation[1];
v[2] = fTranslation[2];
memcpy(fTranslation,v,kN3);
}
void TGeoCombiTrans::ReflectX(Bool_t leftside, Bool_t rotonly)
{
if (leftside && !rotonly) fTranslation[0] = -fTranslation[0];
if (!fRotation || !TestBit(kGeoMatrixOwned)) {
if (fRotation) fRotation = new TGeoRotation(*fRotation);
else fRotation = new TGeoRotation();
SetBit(kGeoMatrixOwned);
}
SetBit(kGeoRotation);
fRotation->ReflectX(leftside);
SetBit(kGeoReflection, !IsReflection());
}
void TGeoCombiTrans::ReflectY(Bool_t leftside, Bool_t rotonly)
{
if (leftside && !rotonly) fTranslation[1] = -fTranslation[1];
if (!fRotation || !TestBit(kGeoMatrixOwned)) {
if (fRotation) fRotation = new TGeoRotation(*fRotation);
else fRotation = new TGeoRotation();
SetBit(kGeoMatrixOwned);
}
SetBit(kGeoRotation);
fRotation->ReflectY(leftside);
SetBit(kGeoReflection, !IsReflection());
}
void TGeoCombiTrans::ReflectZ(Bool_t leftside, Bool_t rotonly)
{
if (leftside && !rotonly) fTranslation[2] = -fTranslation[2];
if (!fRotation || !TestBit(kGeoMatrixOwned)) {
if (fRotation) fRotation = new TGeoRotation(*fRotation);
else fRotation = new TGeoRotation();
SetBit(kGeoMatrixOwned);
}
SetBit(kGeoRotation);
fRotation->ReflectZ(leftside);
SetBit(kGeoReflection, !IsReflection());
}
void TGeoCombiTrans::SavePrimitive(ostream &out, Option_t *option )
{
if (TestBit(kGeoSavePrimitive)) return;
out << " // Combi transformation: " << GetName() << endl;
out << " dx = " << fTranslation[0] << ";" << endl;
out << " dy = " << fTranslation[1] << ";" << endl;
out << " dz = " << fTranslation[2] << ";" << endl;
if (fRotation && fRotation->IsRotation()) {
fRotation->SavePrimitive(out,option);
out << " " << GetPointerName() << " = new TGeoCombiTrans(\"" << GetName() << "\", dx,dy,dz,";
out << fRotation->GetPointerName() << ");" << endl;
} else {
out << " " << GetPointerName() << " = new TGeoCombiTrans(\"" << GetName() << "\");" << endl;
out << " " << GetPointerName() << "->SetTranslation(dx,dy,dz);" << endl;
}
TObject::SetBit(kGeoSavePrimitive);
}
void TGeoCombiTrans::SetRotation(const TGeoRotation *rot)
{
if (fRotation && TestBit(kGeoMatrixOwned)) delete fRotation;
fRotation = 0;
ResetBit(TGeoMatrix::kGeoMatrixOwned);
ResetBit(kGeoRotation);
ResetBit(kGeoReflection);
if (!rot) return;
if (!rot->IsRotation()) return;
SetBit(kGeoRotation);
SetBit(kGeoReflection, rot->TestBit(kGeoReflection));
TGeoRotation *rr = (TGeoRotation*)rot;
fRotation = rr;
}
void TGeoCombiTrans::SetRotation(const TGeoRotation &rot)
{
if (fRotation && TestBit(kGeoMatrixOwned)) delete fRotation;
fRotation = 0;
if (!rot.IsRotation()) {
ResetBit(kGeoRotation);
ResetBit(kGeoReflection);
ResetBit(TGeoMatrix::kGeoMatrixOwned);
return;
}
SetBit(kGeoRotation);
SetBit(kGeoReflection, rot.TestBit(kGeoReflection));
fRotation = new TGeoRotation(rot);
SetBit(kGeoMatrixOwned);
}
void TGeoCombiTrans::SetTranslation(const TGeoTranslation &tr)
{
if (tr.IsTranslation()) {
SetBit(kGeoTranslation);
const Double_t *trans = tr.GetTranslation();
memcpy(fTranslation, trans, kN3);
} else {
if (!IsTranslation()) return;
memset(fTranslation, 0, kN3);
ResetBit(kGeoTranslation);
}
}
void TGeoCombiTrans::SetTranslation(Double_t dx, Double_t dy, Double_t dz)
{
fTranslation[0] = dx;
fTranslation[1] = dy;
fTranslation[2] = dz;
if (fTranslation[0] || fTranslation[1] || fTranslation[2]) SetBit(kGeoTranslation);
else ResetBit(kGeoTranslation);
}
void TGeoCombiTrans::SetTranslation(Double_t *vect)
{
fTranslation[0] = vect[0];
fTranslation[1] = vect[1];
fTranslation[2] = vect[2];
if (fTranslation[0] || fTranslation[1] || fTranslation[2]) SetBit(kGeoTranslation);
else ResetBit(kGeoTranslation);
}
const Double_t *TGeoCombiTrans::GetRotationMatrix() const
{
if (!fRotation) return kIdentityMatrix;
return fRotation->GetRotationMatrix();
}
ClassImp(TGeoGenTrans)
TGeoGenTrans::TGeoGenTrans()
{
SetBit(kGeoGenTrans);
for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
for (Int_t j=0; j<3; j++) fScale[j] = 1.0;
fRotation = 0;
}
TGeoGenTrans::TGeoGenTrans(const char *name)
:TGeoCombiTrans(name)
{
SetBit(kGeoGenTrans);
for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
for (Int_t j=0; j<3; j++) fScale[j] = 1.0;
fRotation = 0;
}
TGeoGenTrans::TGeoGenTrans(Double_t dx, Double_t dy, Double_t dz,
Double_t sx, Double_t sy, Double_t sz, TGeoRotation *rot)
:TGeoCombiTrans("")
{
SetBit(kGeoGenTrans);
SetTranslation(dx, dy, dz);
SetScale(sx, sy, sz);
SetRotation(rot);
}
TGeoGenTrans::TGeoGenTrans(const char *name, Double_t dx, Double_t dy, Double_t dz,
Double_t sx, Double_t sy, Double_t sz, TGeoRotation *rot)
:TGeoCombiTrans(name)
{
SetBit(kGeoGenTrans);
SetTranslation(dx, dy, dz);
SetScale(sx, sy, sz);
SetRotation(rot);
}
TGeoGenTrans::~TGeoGenTrans()
{
}
void TGeoGenTrans::Clear(Option_t *)
{
memset(&fTranslation[0], 0, kN3);
memset(&fScale[0], 0, kN3);
if (fRotation) fRotation->Clear();
}
void TGeoGenTrans::SetScale(Double_t sx, Double_t sy, Double_t sz)
{
fScale[0] = sx;
fScale[1] = sy;
fScale[2] = sz;
if (!(Normalize())) {
Error("ctor", "Invalid scale");
return;
}
}
TGeoMatrix& TGeoGenTrans::Inverse() const
{
Error("Inverse", "not implemented");
static TGeoHMatrix h;
h = *this;
return h;
}
Bool_t TGeoGenTrans::Normalize()
{
Double_t normfactor = fScale[0]*fScale[1]*fScale[2];
if (normfactor <= 1E-5) return kFALSE;
for (Int_t i=0; i<3; i++)
fScale[i] /= normfactor;
return kTRUE;
}
ClassImp(TGeoIdentity)
TGeoIdentity::TGeoIdentity()
{
if (!gGeoIdentity) gGeoIdentity = this;
RegisterYourself();
}
TGeoIdentity::TGeoIdentity(const char *name)
:TGeoMatrix(name)
{
if (!gGeoIdentity) gGeoIdentity = this;
RegisterYourself();
}
TGeoMatrix &TGeoIdentity::Inverse() const
{
return *gGeoIdentity;
}
ClassImp(TGeoHMatrix)
TGeoHMatrix::TGeoHMatrix()
{
memset(&fTranslation[0], 0, kN3);
memcpy(fRotationMatrix,kIdentityMatrix,kN9);
memcpy(fScale,kUnitScale,kN3);
}
TGeoHMatrix::TGeoHMatrix(const char* name)
:TGeoMatrix(name)
{
memset(&fTranslation[0], 0, kN3);
memcpy(fRotationMatrix,kIdentityMatrix,kN9);
memcpy(fScale,kUnitScale,kN3);
}
TGeoHMatrix::TGeoHMatrix(const TGeoMatrix &matrix)
:TGeoMatrix(matrix)
{
if (matrix.IsTranslation()) {
SetBit(kGeoTranslation);
SetTranslation(matrix.GetTranslation());
} else {
memset(&fTranslation[0], 0, kN3);
}
if (matrix.IsRotation()) {
SetBit(kGeoRotation);
memcpy(fRotationMatrix,matrix.GetRotationMatrix(),kN9);
} else {
memcpy(fRotationMatrix,kIdentityMatrix,kN9);
}
if (matrix.IsScale()) {
SetBit(kGeoScale);
memcpy(fScale,matrix.GetScale(),kN3);
} else {
memcpy(fScale,kUnitScale,kN3);
}
}
TGeoHMatrix::~TGeoHMatrix()
{
}
TGeoHMatrix &TGeoHMatrix::operator=(const TGeoMatrix *matrix)
{
if (matrix == this) return *this;
Clear();
if (matrix == 0) return *this;
TGeoMatrix::operator=(*matrix);
if (matrix->IsIdentity()) return *this;
if (matrix->IsTranslation()) {
SetBit(kGeoTranslation);
memcpy(fTranslation,matrix->GetTranslation(),kN3);
}
if (matrix->IsRotation()) {
SetBit(kGeoRotation);
memcpy(fRotationMatrix,matrix->GetRotationMatrix(),kN9);
}
if (matrix->IsScale()) {
SetBit(kGeoScale);
memcpy(fScale,matrix->GetScale(),kN3);
}
return *this;
}
TGeoHMatrix &TGeoHMatrix::operator=(const TGeoMatrix &matrix)
{
if (&matrix == this) return *this;
Clear();
TGeoMatrix::operator=(matrix);
if (matrix.IsIdentity()) return *this;
if (matrix.IsTranslation()) {
SetBit(kGeoTranslation);
memcpy(fTranslation,matrix.GetTranslation(),kN3);
} else {
memcpy(fTranslation,kNullVector,kN3);
}
if (matrix.IsRotation()) {
SetBit(kGeoRotation);
memcpy(fRotationMatrix,matrix.GetRotationMatrix(),kN9);
} else {
memcpy(fRotationMatrix,kIdentityMatrix,kN9);
}
if (matrix.IsScale()) {
SetBit(kGeoScale);
memcpy(fScale,matrix.GetScale(),kN3);
} else {
memcpy(fScale,kUnitScale,kN3);
}
return *this;
}
void TGeoHMatrix::Clear(Option_t *)
{
SetBit(kGeoReflection, kFALSE);
if (IsIdentity()) return;
if (IsTranslation()) {
ResetBit(kGeoTranslation);
memcpy(fTranslation,kNullVector,kN3);
}
if (IsRotation()) {
ResetBit(kGeoRotation);
memcpy(fRotationMatrix,kIdentityMatrix,kN9);
}
if (IsScale()) {
ResetBit(kGeoScale);
memcpy(fScale,kUnitScale,kN3);
}
}
TGeoMatrix *TGeoHMatrix::MakeClone() const
{
TGeoMatrix *matrix = new TGeoHMatrix(*this);
return matrix;
}
TGeoMatrix& TGeoHMatrix::Inverse() const
{
static TGeoHMatrix h;
h = *this;
if (IsTranslation()) {
Double_t tr[3];
tr[0] = -fTranslation[0]*fRotationMatrix[0] - fTranslation[1]*fRotationMatrix[3] - fTranslation[2]*fRotationMatrix[6];
tr[1] = -fTranslation[0]*fRotationMatrix[1] - fTranslation[1]*fRotationMatrix[4] - fTranslation[2]*fRotationMatrix[7];
tr[2] = -fTranslation[0]*fRotationMatrix[2] - fTranslation[1]*fRotationMatrix[5] - fTranslation[2]*fRotationMatrix[8];
h.SetTranslation(tr);
}
if (IsRotation()) {
Double_t newrot[9];
newrot[0] = fRotationMatrix[0];
newrot[1] = fRotationMatrix[3];
newrot[2] = fRotationMatrix[6];
newrot[3] = fRotationMatrix[1];
newrot[4] = fRotationMatrix[4];
newrot[5] = fRotationMatrix[7];
newrot[6] = fRotationMatrix[2];
newrot[7] = fRotationMatrix[5];
newrot[8] = fRotationMatrix[8];
h.SetRotation(newrot);
}
if (IsScale()) {
Double_t sc[3];
sc[0] = 1./fScale[0];
sc[1] = 1./fScale[1];
sc[2] = 1./fScale[2];
h.SetScale(sc);
}
return h;
}
Double_t TGeoHMatrix::Determinant() const
{
Double_t
det = fRotationMatrix[0]*fRotationMatrix[4]*fRotationMatrix[8] +
fRotationMatrix[3]*fRotationMatrix[7]*fRotationMatrix[2] +
fRotationMatrix[6]*fRotationMatrix[1]*fRotationMatrix[5] -
fRotationMatrix[2]*fRotationMatrix[4]*fRotationMatrix[6] -
fRotationMatrix[5]*fRotationMatrix[7]*fRotationMatrix[0] -
fRotationMatrix[8]*fRotationMatrix[1]*fRotationMatrix[3];
return det;
}
void TGeoHMatrix::Multiply(const TGeoMatrix *right)
{
if (right->IsIdentity()) return;
const Double_t *r_tra = right->GetTranslation();
const Double_t *r_rot = right->GetRotationMatrix();
const Double_t *r_scl = right->GetScale();
if (IsIdentity()) {
if (right->IsRotation()) {
SetBit(kGeoRotation);
memcpy(fRotationMatrix,r_rot,kN9);
if (right->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
}
if (right->IsScale()) {
SetBit(kGeoScale);
memcpy(fScale,r_scl,kN3);
}
if (right->IsTranslation()) {
SetBit(kGeoTranslation);
memcpy(fTranslation,r_tra,kN3);
}
return;
}
Int_t i, j;
Double_t new_rot[9];
if (right->IsRotation()) {
SetBit(kGeoRotation);
if (right->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
}
if (right->IsScale()) SetBit(kGeoScale);
if (right->IsTranslation()) SetBit(kGeoTranslation);
if (IsTranslation()) {
for (i=0; i<3; i++) {
fTranslation[i] += fRotationMatrix[3*i]*r_tra[0]
+ fRotationMatrix[3*i+1]*r_tra[1]
+ fRotationMatrix[3*i+2]*r_tra[2];
}
}
if (IsRotation()) {
for (i=0; i<3; i++) {
for (j=0; j<3; j++) {
new_rot[3*i+j] = fRotationMatrix[3*i]*r_rot[j] +
fRotationMatrix[3*i+1]*r_rot[3+j] +
fRotationMatrix[3*i+2]*r_rot[6+j];
}
}
memcpy(fRotationMatrix,new_rot,kN9);
}
if (IsScale()) {
for (i=0; i<3; i++) fScale[i] *= r_scl[i];
}
}
void TGeoHMatrix::MultiplyLeft(const TGeoMatrix *left)
{
if (left == gGeoIdentity) return;
const Double_t *l_tra = left->GetTranslation();
const Double_t *l_rot = left->GetRotationMatrix();
const Double_t *l_scl = left->GetScale();
if (IsIdentity()) {
if (left->IsRotation()) {
if (left->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
SetBit(kGeoRotation);
memcpy(fRotationMatrix,l_rot,kN9);
}
if (left->IsScale()) {
SetBit(kGeoScale);
memcpy(fScale,l_scl,kN3);
}
if (left->IsTranslation()) {
SetBit(kGeoTranslation);
memcpy(fTranslation,l_tra,kN3);
}
return;
}
Int_t i, j;
Double_t new_tra[3];
Double_t new_rot[9];
if (left->IsRotation()) {
SetBit(kGeoRotation);
if (left->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
}
if (left->IsScale()) SetBit(kGeoScale);
if (left->IsTranslation()) SetBit(kGeoTranslation);
if (IsTranslation()) {
for (i=0; i<3; i++) {
new_tra[i] = l_tra[i]
+ l_rot[3*i]* fTranslation[0]
+ l_rot[3*i+1]*fTranslation[1]
+ l_rot[3*i+2]*fTranslation[2];
}
memcpy(fTranslation,new_tra,kN3);
}
if (IsRotation()) {
for (i=0; i<3; i++) {
for (j=0; j<3; j++) {
new_rot[3*i+j] = l_rot[3*i]*fRotationMatrix[j] +
l_rot[3*i+1]*fRotationMatrix[3+j] +
l_rot[3*i+2]*fRotationMatrix[6+j];
}
}
memcpy(fRotationMatrix,new_rot,kN9);
}
if (IsScale()) {
for (i=0; i<3; i++) fScale[i] *= l_scl[i];
}
}
void TGeoHMatrix::RotateX(Double_t angle)
{
SetBit(kGeoRotation);
Double_t phi = angle*TMath::DegToRad();
Double_t c = TMath::Cos(phi);
Double_t s = TMath::Sin(phi);
Double_t v[9];
v[0] = fRotationMatrix[0];
v[1] = fRotationMatrix[1];
v[2] = fRotationMatrix[2];
v[3] = c*fRotationMatrix[3]-s*fRotationMatrix[6];
v[4] = c*fRotationMatrix[4]-s*fRotationMatrix[7];
v[5] = c*fRotationMatrix[5]-s*fRotationMatrix[8];
v[6] = s*fRotationMatrix[3]+c*fRotationMatrix[6];
v[7] = s*fRotationMatrix[4]+c*fRotationMatrix[7];
v[8] = s*fRotationMatrix[5]+c*fRotationMatrix[8];
memcpy(fRotationMatrix, v, kN9);
v[0] = fTranslation[0];
v[1] = c*fTranslation[1]-s*fTranslation[2];
v[2] = s*fTranslation[1]+c*fTranslation[2];
memcpy(fTranslation,v,kN3);
}
void TGeoHMatrix::RotateY(Double_t angle)
{
SetBit(kGeoRotation);
Double_t phi = angle*TMath::DegToRad();
Double_t c = TMath::Cos(phi);
Double_t s = TMath::Sin(phi);
Double_t v[9];
v[0] = c*fRotationMatrix[0]+s*fRotationMatrix[6];
v[1] = c*fRotationMatrix[1]+s*fRotationMatrix[7];
v[2] = c*fRotationMatrix[2]+s*fRotationMatrix[8];
v[3] = fRotationMatrix[3];
v[4] = fRotationMatrix[4];
v[5] = fRotationMatrix[5];
v[6] = -s*fRotationMatrix[0]+c*fRotationMatrix[6];
v[7] = -s*fRotationMatrix[1]+c*fRotationMatrix[7];
v[8] = -s*fRotationMatrix[2]+c*fRotationMatrix[8];
memcpy(fRotationMatrix, v, kN9);
v[0] = c*fTranslation[0]+s*fTranslation[2];
v[1] = fTranslation[1];
v[2] = -s*fTranslation[0]+c*fTranslation[2];
memcpy(fTranslation,v,kN3);
}
void TGeoHMatrix::RotateZ(Double_t angle)
{
SetBit(kGeoRotation);
Double_t phi = angle*TMath::DegToRad();
Double_t c = TMath::Cos(phi);
Double_t s = TMath::Sin(phi);
Double_t v[9];
v[0] = c*fRotationMatrix[0]-s*fRotationMatrix[3];
v[1] = c*fRotationMatrix[1]-s*fRotationMatrix[4];
v[2] = c*fRotationMatrix[2]-s*fRotationMatrix[5];
v[3] = s*fRotationMatrix[0]+c*fRotationMatrix[3];
v[4] = s*fRotationMatrix[1]+c*fRotationMatrix[4];
v[5] = s*fRotationMatrix[2]+c*fRotationMatrix[5];
v[6] = fRotationMatrix[6];
v[7] = fRotationMatrix[7];
v[8] = fRotationMatrix[8];
memcpy(&fRotationMatrix[0],v,kN9);
v[0] = c*fTranslation[0]-s*fTranslation[1];
v[1] = s*fTranslation[0]+c*fTranslation[1];
v[2] = fTranslation[2];
memcpy(fTranslation,v,kN3);
}
void TGeoHMatrix::ReflectX(Bool_t leftside, Bool_t rotonly)
{
if (leftside && !rotonly) fTranslation[0] = -fTranslation[0];
if (leftside) {
fRotationMatrix[0]=-fRotationMatrix[0];
fRotationMatrix[1]=-fRotationMatrix[1];
fRotationMatrix[2]=-fRotationMatrix[2];
} else {
fRotationMatrix[0]=-fRotationMatrix[0];
fRotationMatrix[3]=-fRotationMatrix[3];
fRotationMatrix[6]=-fRotationMatrix[6];
}
SetBit(kGeoRotation);
SetBit(kGeoReflection, !IsReflection());
}
void TGeoHMatrix::ReflectY(Bool_t leftside, Bool_t rotonly)
{
if (leftside && !rotonly) fTranslation[1] = -fTranslation[1];
if (leftside) {
fRotationMatrix[3]=-fRotationMatrix[3];
fRotationMatrix[4]=-fRotationMatrix[4];
fRotationMatrix[5]=-fRotationMatrix[5];
} else {
fRotationMatrix[1]=-fRotationMatrix[1];
fRotationMatrix[4]=-fRotationMatrix[4];
fRotationMatrix[7]=-fRotationMatrix[7];
}
SetBit(kGeoRotation);
SetBit(kGeoReflection, !IsReflection());
}
void TGeoHMatrix::ReflectZ(Bool_t leftside, Bool_t rotonly)
{
if (leftside && !rotonly) fTranslation[2] = -fTranslation[2];
if (leftside) {
fRotationMatrix[6]=-fRotationMatrix[6];
fRotationMatrix[7]=-fRotationMatrix[7];
fRotationMatrix[8]=-fRotationMatrix[8];
} else {
fRotationMatrix[2]=-fRotationMatrix[2];
fRotationMatrix[5]=-fRotationMatrix[5];
fRotationMatrix[8]=-fRotationMatrix[8];
}
SetBit(kGeoRotation);
SetBit(kGeoReflection, !IsReflection());
}
void TGeoHMatrix::SavePrimitive(ostream &out, Option_t * )
{
if (TestBit(kGeoSavePrimitive)) return;
const Double_t *tr = fTranslation;
const Double_t *rot = fRotationMatrix;
out << " // HMatrix: " << GetName() << endl;
out << " tr[0] = " << tr[0] << "; " << "tr[1] = " << tr[1] << "; " << "tr[2] = " << tr[2] << ";" << endl;
out << " rot[0] =" << rot[0] << "; " << "rot[1] = " << rot[1] << "; " << "rot[2] = " << rot[2] << ";" << endl;
out << " rot[3] =" << rot[3] << "; " << "rot[4] = " << rot[4] << "; " << "rot[5] = " << rot[5] << ";" << endl;
out << " rot[6] =" << rot[6] << "; " << "rot[7] = " << rot[7] << "; " << "rot[8] = " << rot[8] << ";" << endl;
char *name = GetPointerName();
out << " TGeoHMatrix *" << name << " = new TGeoHMatrix(\"" << GetName() << "\");" << endl;
out << " " << name << "->SetTranslation(tr);" << endl;
out << " " << name << "->SetRotation(rot);" << endl;
if (IsTranslation()) out << " " << name << "->SetBit(TGeoMatrix::kGeoTranslation);" << endl;
if (IsRotation()) out << " " << name << "->SetBit(TGeoMatrix::kGeoRotation);" << endl;
if (IsReflection()) out << " " << name << "->SetBit(TGeoMatrix::kGeoReflection);" << endl;
TObject::SetBit(kGeoSavePrimitive);
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.