STEER/AliAlignObj.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-----------------------------------------------------------------
  19 //  Implementation of the alignment object class, holding the alignment
  20 //  constants for a single volume, through the abstract class AliAlignObj.
  21 //  From it two derived concrete representation of alignment object class
  22 //  (AliAlignObjParams, AliAlignObjMatrix) are derived in separate files.
  23 //-----------------------------------------------------------------
  24
  25 #include <TGeoManager.h>
  26 #include <TGeoMatrix.h>
  27 #include <TGeoPhysicalNode.h>
  28 #include <TGeoOverlap.h>
  29 #include <TMath.h>
  30
  31 #include "AliAlignObj.h"
  32 #include "AliTrackPointArray.h"
  33 #include "AliLog.h"
  34
  35 ClassImp(AliAlignObj)
  36
  37 //_____________________________________________________________________________
  38 AliAlignObj::AliAlignObj():
  39   TObject(),
  40   fVolPath(),
  41   fVolUID(0)
  42 {
  43   // default constructor
  44   for(Int_t i=0; i<6; i++) fDiag[i]=-999.;
  45   for(Int_t i=0; i<15; i++) fODia[i]=-999.;
  46 }
  47
  48 //_____________________________________________________________________________
  49 AliAlignObj::AliAlignObj(const char* symname, UShort_t voluid) :
  50   TObject(),
  51   fVolPath(symname),
  52   fVolUID(voluid)
  53 {
  54   // standard constructor
  55   //
  56   for(Int_t i=0; i<6; i++) fDiag[i]=-999.;
  57   for(Int_t i=0; i<15; i++) fODia[i]=-999.;
  58 }
  59
  60 //_____________________________________________________________________________
  61 AliAlignObj::AliAlignObj(const char* symname, UShort_t voluid, Double_t* cmat) :
  62   TObject(),
  63   fVolPath(symname),
  64   fVolUID(voluid)
  65 {
  66   // standard constructor
  67   //
  68   SetCorrMatrix(cmat);
  69 }
  70
  71 //_____________________________________________________________________________
  72 AliAlignObj::AliAlignObj(const AliAlignObj& theAlignObj) :
  73   TObject(theAlignObj),
  74   fVolPath(theAlignObj.GetSymName()),
  75   fVolUID(theAlignObj.GetVolUID())
  76 {
  77   //copy constructor
  78   for(Int_t i=0; i<6; i++) fDiag[i]=theAlignObj.fDiag[i];
  79   for(Int_t i=0; i<15; i++) fODia[i]=theAlignObj.fODia[i];
  80 }
  81
  82 //_____________________________________________________________________________
  83 AliAlignObj &AliAlignObj::operator =(const AliAlignObj& theAlignObj)
  84 {
  85   // assignment operator
  86   if(this==&theAlignObj) return *this;
  87   fVolPath = theAlignObj.GetSymName();
  88   fVolUID = theAlignObj.GetVolUID();
  89   for(Int_t i=0; i<6; i++) fDiag[i]=theAlignObj.fDiag[i];
  90   for(Int_t i=0; i<15; i++) fODia[i]=theAlignObj.fODia[i];
  91   return *this;
  92 }
  93
  94 //_____________________________________________________________________________
  95 AliAlignObj &AliAlignObj::operator*=(const AliAlignObj& theAlignObj)
  96 {
  97   // multiplication operator
  98   // The operator can be used to 'combine'
  99   // two alignment objects
 100   TGeoHMatrix m1;
 101   GetMatrix(m1);
 102   TGeoHMatrix m2;
 103   theAlignObj.GetMatrix(m2);
 104   m1.MultiplyLeft(&m2);
 105   SetMatrix(m1);
 106   // temporary solution: the covariance matrix of the resulting combined object
 107   // is set equal to the covariance matrix of the right operand
 108   // (not to be used for combining alignment objects for different levels)
 109   for(Int_t i=0; i<6; i++)  fDiag[i] = theAlignObj.fDiag[i];
 110   for(Int_t i=0; i<15; i++)  fODia[i] = theAlignObj.fODia[i];
 111   return *this;
 112 }
 113
 114 //_____________________________________________________________________________
 115 AliAlignObj::~AliAlignObj()
 116 {
 117   // dummy destructor
 118 }
 119
 120 //_____________________________________________________________________________
 121 void AliAlignObj::SetVolUID(AliGeomManager::ELayerID detId, Int_t modId)
 122 {
 123   // From detector name and module number (according to detector numbering)
 124   // build fVolUID, unique numerical identity of that volume inside ALICE
 125   // fVolUID is 16 bits, first 5 reserved for detID (32 possible values),
 126   // remaining 11 for module ID inside det (2048 possible values).
 127   //
 128   fVolUID = AliGeomManager::LayerToVolUID(detId,modId);
 129 }
 130
 131 //_____________________________________________________________________________
 132 void AliAlignObj::GetVolUID(AliGeomManager::ELayerID &layerId, Int_t &modId) const
 133 {
 134   // From the fVolUID, unique numerical identity of that volume inside ALICE,
 135   // (voluid is 16 bits, first 5 reserved for layerID (32 possible values),
 136   // remaining 11 for module ID inside det (2048 possible values)), sets
 137   // the argument layerId to the identity of the layer to which that volume
 138   // belongs and sets the argument modId to the identity of that volume
 139   // internally to the layer.
 140   //
 141   layerId = AliGeomManager::VolUIDToLayer(fVolUID,modId);
 142 }
 143
 144 //_____________________________________________________________________________
 145 Bool_t AliAlignObj::GetPars(Double_t tr[], Double_t angles[]) const
 146 {
 147   GetTranslation(tr);
 148   return GetAngles(angles);
 149 }
 150
 151 //_____________________________________________________________________________
 152 Int_t AliAlignObj::GetLevel() const
 153 {
 154   // Return the geometry level of the alignable volume to which
 155   // the alignment object is associated; this is the number of
 156   // slashes in the corresponding volume path
 157   //
 158   if(!gGeoManager){
 159     AliWarning("gGeoManager doesn't exist or it is still open: unable to return meaningful level value.");
 160     return (-1);
 161   }
 162   const char* symname = GetSymName();
 163   const char* path;
 164   TGeoPNEntry* pne = gGeoManager->GetAlignableEntry(symname);
 165   if(pne){
 166     path = pne->GetTitle();
 167   }else{
 168     path = symname;
 169   }
 170
 171   TString pathStr = path;
 172   if(pathStr[0]!='/') pathStr.Prepend('/');
 173   return pathStr.CountChar('/');
 174 }
 175
 176 //_____________________________________________________________________________
 177 Int_t AliAlignObj::Compare(const TObject *obj) const
 178 {
 179   // Compare the levels of two
 180   // alignment objects
 181   // Used in the sorting during
 182   // the application of alignment
 183   // objects to the geometry
 184   //
 185   Int_t level = GetLevel();
 186   Int_t level2 = ((AliAlignObj *)obj)->GetLevel();
 187   if (level == level2)
 188     return 0;
 189   else
 190     return ((level > level2) ? 1 : -1);
 191 }
 192
 193 //______________________________________________________________________________
 194 void AliAlignObj::GetCovMatrix(Double_t *cmat) const
 195 {
 196   // Fills the cmat argument with the coefficients of the external cov matrix (21 elements)
 197   // calculating them from the correlation matrix data member
 198   //
 199
 200   for(Int_t i=0; i<6; ++i) {
 201     // Off diagonal elements
 202     for(Int_t j=0; j<i; ++j) {
 203       cmat[i*(i+1)/2+j] = (fDiag[j] >= 0. && fDiag[i] >= 0.) ? fODia[(i-1)*i/2+j]*fDiag[j]*fDiag[i]: -999.;
 204     }
 205
 206     // Diagonal elements
 207     cmat[i*(i+1)/2+i] = (fDiag[i] >= 0.) ? fDiag[i]*fDiag[i] : -999.;
 208   }
 209
 210   return;
 211 }
 212
 213 //______________________________________________________________________________
 214 void AliAlignObj::GetCovMatrix(TMatrixDSym& mcov) const
 215 {
 216   // Fills the matrix m passed as argument as the covariance matrix calculated
 217   // from the coefficients of the reduced covariance matrix data members
 218   //
 219
 220   for(Int_t i=0; i<6; ++i) {
 221     // Off diagonal elements
 222     for(Int_t j=0; j<i; ++j) {
 223       mcov(j,i) = mcov(i,j) = (fDiag[j] >= 0. && fDiag[i] >= 0.) ? fODia[(i-1)*i/2+j]*fDiag[j]*fDiag[i]: -999.;
 224     }
 225
 226     // Diagonal elements
 227     mcov(i,i) = (fDiag[i] >= 0.) ? fDiag[i]*fDiag[i] : -999.;
 228   }
 229
 230 }
 231
 232 //______________________________________________________________________________
 233 Bool_t AliAlignObj::GetLocalCovMatrix(TMatrixDSym& lCov) const
 234 {
 235   // Calculates the covariance matrix (6x6) associated to the six parameters
 236   // defining the current alignment in the global coordinates system (and sets
 237   // in the internal data members) from the covariance matrix (6x6) for the six
 238   // parameters defining the alignment transformation in the local coordinates
 239   // system, passed as an argument.
 240   //
 241   TMatrixD mJ(6,6);// the jacobian of the transformation from local to global parameters
 242   if(!GetJacobian(mJ)) return kFALSE;
 243
 244   TMatrixDSym gCov(6);
 245   GetCovMatrix(gCov);
 246
 247   // Compute the local covariance matrix lcov = mJ^T gcov mJ
 248   TMatrixD gcovJ(gCov,TMatrixD::kMult,mJ);
 249   TMatrixD lCovM(mJ,TMatrixD::kTransposeMult,gcovJ);
 250   // To be done: somehow check that lCovM is close enough to be symmetric
 251   for(Int_t i=0; i<6; i++)
 252   {
 253     lCov(i,i) = lCovM(i,i);
 254     for(Int_t j=i+1; j<6; j++)
 255     {
 256       lCov(i,j)=lCovM(i,j);
 257       lCov(j,i)=lCovM(i,j);
 258     }
 259   }
 260
 261   return kTRUE;
 262
 263 }
 264
 265 //______________________________________________________________________________
 266 Bool_t AliAlignObj::GetLocalCovMatrix(Double_t *lCov) const
 267 {
 268   // Calculates the covariance matrix (6x6) associated to the six parameters
 269   // defining the current alignment in the global coordinates system (and sets
 270   // in the internal data members) from the covariance matrix (6x6) for the six
 271   // parameters defining the alignment transformation in the local coordinates
 272   // system, passed as an argument.
 273   //
 274   TMatrixDSym lCovMatrix(6);
 275   GetLocalCovMatrix(lCovMatrix);
 276
 277   Int_t k=0;
 278   for(Int_t i=0; i<6; i++)
 279     for(Int_t j=i; j<6; j++)
 280     {
 281        lCov[k++] = lCovMatrix(i,j);
 282     }
 283
 284   return kTRUE;
 285 }
 286
 287 //______________________________________________________________________________
 288 Bool_t AliAlignObj::GetJacobian(TMatrixD& mJ) const
 289 {
 290   // Compute the jacobian J of the transformation of the six local to the six global delta parameters
 291   //
 292   // R00 R01 R02 | (R01Rk2 - R02Rk1)Tk  (R02Rk0 - R00Rk2)Tk  (R00Rk1 - R01Rk0)Tk
 293   // R00 R01 R02 | (R11Rk2 - R12Rk1)Tk  (R12Rk0 - R10Rk2)Tk  (R10Rk1 - R11Rk0)Tk
 294   // R00 R01 R02 | (R21Rk2 - R22Rk1)Tk  (R22Rk0 - R20Rk2)Tk  (R20Rk1 - R21Rk0)Tk
 295   //  -  -   -   -   -   -   -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
 296   //  0   0   0  |   R11R22 - R12R21      R12R20 - R10R22      R10R21 - R11R20
 297   //  0   0   0  |   R21R02 - R22R01      R22R00 - R20R02      R20R01 - R21R00
 298   //  0   0   0  |   R01R12 - R02R11      R02R10 - R00R12      R00R11 - R01R10
 299   //
 300   if (!gGeoManager || !gGeoManager->IsClosed()) {
 301     AliError("Can't compute the global covariance matrix from the local one without an open geometry!");
 302     return kFALSE;
 303   }
 304
 305   const char* symname = GetSymName();
 306   TGeoPhysicalNode* node;
 307   TGeoPNEntry* pne = gGeoManager->GetAlignableEntry(symname);
 308   if(pne){
 309     if(!pne->GetPhysicalNode()){
 310       node = gGeoManager->MakeAlignablePN(pne);
 311     }else{
 312       node = pne->GetPhysicalNode();
 313     }
 314   }else{
 315     AliWarning(Form("The symbolic volume name %s does not correspond to a physical entry. Using it as volume path!",symname));
 316     node = (TGeoPhysicalNode*) gGeoManager->MakePhysicalNode(symname);
 317   }
 318
 319   if (!node) {
 320     AliError(Form("Volume name or path %s not valid!",symname));
 321     return kFALSE;
 322   }
 323
 324   TGeoHMatrix gm; //global matrix
 325   gm = *node->GetMatrix();
 326   Double_t *tr  = gm.GetTranslation();
 327   Double_t *rot = gm.GetRotationMatrix();
 328
 329   TGeoHMatrix m; // global delta transformation matrix
 330   GetMatrix(m);
 331   // We should probably check that it's sufficinetly close to identity
 332   // if it's not return because the "small angles" approximation cannot hold
 333
 334   // 3x3 upper left part (global shifts derived w.r.t. local shifts)
 335   for(Int_t i=0; i<3; i++)
 336   {
 337     for(Int_t j=0; j<3; j++)
 338     {
 339       mJ(i,j) = rot[i+3*j];
 340     }
 341   }
 342
 343   // 3x3 lower left part (global angles derived w.r.t. local shifts)
 344   for(Int_t i=0; i<3; i++)
 345   {
 346     for(Int_t j=0; j<3; j++)
 347     {
 348       mJ(i+3,j) = 0.;
 349     }
 350   }
 351
 352   // 3x3 upper right part (global shifts derived w.r.t. local angles)
 353   for(Int_t i=0; i<3; i++)
 354   {
 355     for(Int_t j=0; j<3; j++)
 356     {
 357       Double_t mEl = 0.;
 358       Int_t b = (j+1)%3;
 359       Int_t d = (j+2)%3;
 360       for(Int_t k=0; k<3; k++)
 361       {
 362         mEl += (rot[3*i+b]*rot[3*k+d])*tr[k]-(rot[3*i+d]*rot[3*k+b])*tr[k];
 363       }
 364       mJ(i,j+3) = mEl;
 365     }
 366   }
 367
 368   // 3x3 lower right part (global angles derived w.r.t. local angles)
 369   for(Int_t i=0; i<3; i++)
 370     for(Int_t j=0; j<3; j++)
 371     {
 372       Int_t a = (i+1)%3;
 373       Int_t b = (j+1)%3;
 374       Int_t c = (i+2)%3;
 375       Int_t d = (j+2)%3;
 376       mJ(i+3,j+3) = rot[3*a+b]*rot[3*c+d]-rot[3*a+d]*rot[3*c+b];
 377     }
 378
 379   return kTRUE;
 380
 381 }
 382
 383 //______________________________________________________________________________
 384 Bool_t AliAlignObj::SetFromLocalCov(TMatrixDSym& lCov)
 385 {
 386   // Calculates the covariance matrix (6x6) associated to the six parameters
 387   // defining the current alignment in the global coordinates system (and sets
 388   // in the internal data members) from the covariance matrix (6x6) for the six
 389   // parameters defining the alignment transformation in the local coordinates
 390   // system, passed as an argument.
 391   //
 392   TMatrixD mJ(6,6);// the jacobian of the transformation from local to global parameters
 393   if(!GetJacobian(mJ)) return kFALSE;
 394
 395   // Compute the global covariance matrix gcov = mJ lcov mJ'
 396   TMatrixD trJ(TMatrixD::kTransposed, mJ);
 397   TMatrixD lcovTrJ(lCov,TMatrixD::kMult,trJ);
 398   TMatrixD gCovM(mJ,TMatrixD::kMult,lcovTrJ);
 399   // To be done: somehow check that gCovM is close enough to be symmetric
 400   TMatrixDSym gCov(6);
 401   for(Int_t i=0; i<6; i++)
 402   {
 403     gCov(i,i) = gCovM(i,i);
 404     for(Int_t j=i+1; j<6; j++)
 405     {
 406       gCov(i,j)=gCovM(i,j);
 407       gCov(j,i)=gCovM(i,j);
 408     }
 409   }
 410   SetCorrMatrix(gCov);
 411
 412   return kTRUE;
 413
 414 }
 415
 416 //______________________________________________________________________________
 417 Bool_t AliAlignObj::SetFromLocalCov(Double_t *lCov)
 418 {
 419   // Calculates the covariance matrix (6x6) associated to the six parameters
 420   // defining the current alignment in the global coordinates system, and sets
 421   // in the internal data members, from the 21 coefficients, passed as argument,
 422   // of the covariance matrix (6x6) for the six parameters defining the
 423   // alignment transformation in the local coordinates system.
 424   //
 425   TMatrixDSym lCovMatrix(6);
 426
 427   Int_t k=0;
 428   for(Int_t i=0; i<6; i++)
 429     for(Int_t j=i; j<6; j++)
 430     {
 431       lCovMatrix(i,j) = lCov[k++];
 432       if(j!=i) lCovMatrix(j,i) = lCovMatrix(i,j);
 433     }
 434
 435   return SetFromLocalCov(lCovMatrix);
 436
 437 }
 438
 439
 440 //______________________________________________________________________________
 441 void AliAlignObj::SetCorrMatrix(Double_t *cmat)
 442 {
 443   // Sets the correlation matrix data member from the coefficients of the external covariance
 444   // matrix (21 elements passed as argument).
 445   //
 446   if(cmat) {
 447
 448     // Diagonal elements first
 449     for(Int_t i=0; i<6; ++i) {
 450       fDiag[i] = (cmat[i*(i+1)/2+i] >= 0.) ? TMath::Sqrt(cmat[i*(i+1)/2+i]) : -999.;
 451     }
 452
 453     // ... then the ones off diagonal
 454     for(Int_t i=0; i<6; ++i)
 455       // Off diagonal elements
 456       for(Int_t j=0; j<i; ++j) {
 457         fODia[(i-1)*i/2+j] = (fDiag[i] > 0. && fDiag[j] > 0.) ? cmat[i*(i+1)/2+j]/(fDiag[j]*fDiag[i]) : 0.;       // check for division by zero (due to diagonal element of 0) and for fDiag != -999. (due to negative input diagonal element).
 458         if (fODia[(i-1)*i/2+j]>1.)  fODia[(i-1)*i/2+j] =  1.; // check upper boundary
 459         if (fODia[(i-1)*i/2+j]<-1.) fODia[(i-1)*i/2+j] = -1.; // check lower boundary
 460       }
 461   } else {
 462     for(Int_t i=0; i< 6; ++i) fDiag[i]=-999.;
 463     for(Int_t i=0; i< 6*(6-1)/2; ++i) fODia[i]=0.;
 464   }
 465
 466   return;
 467 }
 468
 469 //______________________________________________________________________________
 470 void AliAlignObj::SetCorrMatrix(TMatrixDSym& mcov)
 471 {
 472   // Sets the correlation matrix data member from the covariance matrix mcov passed
 473   // passed as argument.
 474   //
 475   if(mcov.IsValid()) {
 476
 477     // Diagonal elements first
 478     for(Int_t i=0; i<6; ++i) {
 479       fDiag[i] = (mcov(i,i) >= 0.) ? TMath::Sqrt(mcov(i,i)) : -999.;
 480     }
 481
 482     // ... then the ones off diagonal
 483     for(Int_t i=0; i<6; ++i)
 484       // Off diagonal elements
 485       for(Int_t j=0; j<i; ++j) {
 486         fODia[(i-1)*i/2+j] = (fDiag[i] > 0. && fDiag[j] > 0.) ? mcov(i,j)/(fDiag[j]*fDiag[i]) : 0.;       // check for division by zero (due to diagonal element of 0) and for fDiag != -999. (due to negative input diagonal element).
 487         if (fODia[(i-1)*i/2+j]>1.)  fODia[(i-1)*i/2+j] =  1.; // check upper boundary
 488         if (fODia[(i-1)*i/2+j]<-1.) fODia[(i-1)*i/2+j] = -1.; // check lower boundary
 489       }
 490   } else {
 491     for(Int_t i=0; i< 6; ++i) fDiag[i]=-999.;
 492     for(Int_t i=0; i< 6*(6-1)/2; ++i) fODia[i]=0.;
 493   }
 494
 495   return;
 496 }
 497
 498 //_____________________________________________________________________________
 499 void AliAlignObj::AnglesToMatrix(const Double_t *angles, Double_t *rot) const
 500 {
 501   // Calculates the rotation matrix using the
 502   // Euler angles in "x y z" notation
 503   //
 504   Double_t degrad = TMath::DegToRad();
 505   Double_t sinpsi = TMath::Sin(degrad*angles[0]);
 506   Double_t cospsi = TMath::Cos(degrad*angles[0]);
 507   Double_t sinthe = TMath::Sin(degrad*angles[1]);
 508   Double_t costhe = TMath::Cos(degrad*angles[1]);
 509   Double_t sinphi = TMath::Sin(degrad*angles[2]);
 510   Double_t cosphi = TMath::Cos(degrad*angles[2]);
 511
 512   rot[0] =  costhe*cosphi;
 513   rot[1] = -costhe*sinphi;
 514   rot[2] =  sinthe;
 515   rot[3] =  sinpsi*sinthe*cosphi + cospsi*sinphi;
 516   rot[4] = -sinpsi*sinthe*sinphi + cospsi*cosphi;
 517   rot[5] = -costhe*sinpsi;
 518   rot[6] = -cospsi*sinthe*cosphi + sinpsi*sinphi;
 519   rot[7] =  cospsi*sinthe*sinphi + sinpsi*cosphi;
 520   rot[8] =  costhe*cospsi;
 521 }
 522
 523 //_____________________________________________________________________________
 524 Bool_t AliAlignObj::MatrixToAngles(const Double_t *rot, Double_t *angles) const
 525 {
 526   // Calculates the Euler angles in "x y z" notation
 527   // using the rotation matrix
 528   // Returns false in case the rotation angles can not be
 529   // extracted from the matrix
 530   //
 531   if(TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) {
 532     AliError("Failed to extract roll-pitch-yall angles!");
 533     return kFALSE;
 534   }
 535   Double_t raddeg = TMath::RadToDeg();
 536   angles[0]=raddeg*TMath::ATan2(-rot[5],rot[8]);
 537   angles[1]=raddeg*TMath::ASin(rot[2]);
 538   angles[2]=raddeg*TMath::ATan2(-rot[1],rot[0]);
 539   return kTRUE;
 540 }
 541
 542 //______________________________________________________________________________
 543 void AliAlignObj::Transform(AliTrackPoint &p, Bool_t copycov) const
 544 {
 545   // The method transforms the space-point coordinates using the
 546   // transformation matrix provided by the AliAlignObj
 547   // In case the copycov flag is set to kTRUE, the covariance matrix
 548   // of the alignment object is copied into the space-point
 549   //
 550   if (fVolUID != p.GetVolumeID())
 551     AliWarning(Form("Alignment object ID is not equal to the space-point ID (%d != %d)",fVolUID,p.GetVolumeID()));
 552
 553   TGeoHMatrix m;
 554   GetMatrix(m);
 555   Double_t *rot = m.GetRotationMatrix();
 556   Double_t *tr  = m.GetTranslation();
 557
 558   Float_t xyzin[3],xyzout[3];
 559   p.GetXYZ(xyzin);
 560   for (Int_t i = 0; i < 3; i++)
 561     xyzout[i] = tr[i]+
 562                 xyzin[0]*rot[3*i]+
 563                 xyzin[1]*rot[3*i+1]+
 564                 xyzin[2]*rot[3*i+2];
 565   p.SetXYZ(xyzout);
 566
 567   if(copycov){
 568     TMatrixDSym covmat(6);
 569     GetCovMatrix(covmat);
 570     p.SetAlignCovMatrix(covmat);
 571   }
 572
 573 }
 574
 575 //_____________________________________________________________________________
 576 void AliAlignObj::Transform(AliTrackPointArray &array) const
 577 {
 578   // This method is used to transform all the track points
 579   // from the input AliTrackPointArray
 580   //
 581   AliTrackPoint p;
 582   for (Int_t i = 0; i < array.GetNPoints(); i++) {
 583     array.GetPoint(p,i);
 584     Transform(p);
 585     array.AddPoint(i,&p);
 586   }
 587 }
 588
 589 //_____________________________________________________________________________
 590 void AliAlignObj::Print(Option_t *) const
 591 {
 592   // Print the contents of the
 593   // alignment object in angles and
 594   // matrix representations
 595   //
 596   Double_t tr[3];
 597   GetTranslation(tr);
 598   Double_t angles[3];
 599   GetAngles(angles);
 600   TGeoHMatrix m;
 601   GetMatrix(m);
 602   const Double_t *rot = m.GetRotationMatrix();
 603
 604   printf("Volume=%s\n",GetSymName());
 605   if (GetVolUID() != 0) {
 606     AliGeomManager::ELayerID layerId;
 607     Int_t modId;
 608     GetVolUID(layerId,modId);
 609     printf("VolumeID=%d LayerID=%d ( %s ) ModuleID=%d\n", GetVolUID(),layerId,AliGeomManager::LayerName(layerId),modId);
 610   }
 611   printf("%12.8f%12.8f%12.8f    Tx = %12.8f    Psi   = %12.8f\n", rot[0], rot[1], rot[2], tr[0], angles[0]);
 612   printf("%12.8f%12.8f%12.8f    Ty = %12.8f    Theta = %12.8f\n", rot[3], rot[4], rot[5], tr[1], angles[1]);
 613   printf("%12.8f%12.8f%12.8f    Tz = %12.8f    Phi   = %12.8f\n", rot[6], rot[7], rot[8], tr[2], angles[2]);
 614
 615 }
 616
 617 //_____________________________________________________________________________
 618 void AliAlignObj::SetPars(Double_t x, Double_t y, Double_t z,
 619                           Double_t psi, Double_t theta, Double_t phi)
 620 {
 621   // Set the global delta transformation by passing 3 angles (expressed in
 622   // degrees) and 3 shifts (in centimeters)
 623   //
 624   SetTranslation(x,y,z);
 625   SetRotation(psi,theta,phi);
 626 }
 627
 628 //_____________________________________________________________________________
 629 Bool_t AliAlignObj::SetLocalPars(Double_t x, Double_t y, Double_t z,
 630                                  Double_t psi, Double_t theta, Double_t phi)
 631 {
 632   // Set the global delta transformation by passing the parameters
 633   // for the local delta transformation (3 shifts and 3 angles).
 634   // In case that the TGeo was not initialized or not closed,
 635   // returns false and the object parameters are not set.
 636   //
 637   TGeoHMatrix m;
 638   Double_t tr[3] = {x, y, z};
 639   m.SetTranslation(tr);
 640   Double_t angles[3] = {psi, theta, phi};
 641   Double_t rot[9];
 642   AnglesToMatrix(angles,rot);
 643   m.SetRotation(rot);
 644
 645   return SetLocalMatrix(m);
 646
 647 }
 648
 649 //_____________________________________________________________________________
 650 Bool_t AliAlignObj::SetLocalTranslation(Double_t x, Double_t y, Double_t z)
 651 {
 652   // Set the global delta transformation by passing the three shifts giving
 653   // the translation in the local reference system of the alignable
 654   // volume (known by TGeo geometry).
 655   // In case that the TGeo was not initialized or not closed,
 656   // returns false and the object parameters are not set.
 657   //
 658   TGeoHMatrix m;
 659   Double_t tr[3] = {x, y, z};
 660   m.SetTranslation(tr);
 661
 662   return SetLocalMatrix(m);
 663
 664 }
 665
 666 //_____________________________________________________________________________
 667 Bool_t AliAlignObj::SetLocalTranslation(const TGeoMatrix& m)
 668 {
 669   // Set the global delta transformation by passing the matrix of
 670   // the local delta transformation and taking its translational part
 671   // In case that the TGeo was not initialized or not closed,
 672   // returns false and the object parameters are not set.
 673   //
 674   const Double_t* tr = m.GetTranslation();
 675   TGeoHMatrix mtr;
 676   mtr.SetTranslation(tr);
 677
 678   return SetLocalMatrix(mtr);
 679
 680 }
 681
 682 //_____________________________________________________________________________
 683 Bool_t AliAlignObj::SetLocalRotation(Double_t psi, Double_t theta, Double_t phi)
 684 {
 685   // Set the global delta transformation by passing the three angles giving
 686   // the rotation in the local reference system of the alignable
 687   // volume (known by TGeo geometry).
 688   // In case that the TGeo was not initialized or not closed,
 689   // returns false and the object parameters are not set.
 690   //
 691   TGeoHMatrix m;
 692   Double_t angles[3] = {psi, theta, phi};
 693   Double_t rot[9];
 694   AnglesToMatrix(angles,rot);
 695   m.SetRotation(rot);
 696
 697   return SetLocalMatrix(m);
 698
 699 }
 700
 701 //_____________________________________________________________________________
 702 Bool_t AliAlignObj::SetLocalRotation(const TGeoMatrix& m)
 703 {
 704   // Set the global delta transformation by passing the matrix of
 705   // the local delta transformation and taking its rotational part
 706   // In case that the TGeo was not initialized or not closed,
 707   // returns false and the object parameters are not set.
 708   //
 709   TGeoHMatrix rotm;
 710   const Double_t* rot = m.GetRotationMatrix();
 711   rotm.SetRotation(rot);
 712
 713   return SetLocalMatrix(rotm);
 714
 715 }
 716
 717 //_____________________________________________________________________________
 718 Bool_t AliAlignObj::SetLocalMatrix(const TGeoMatrix& m)
 719 {
 720   // Set the global delta transformation by passing the TGeo matrix
 721   // for the local delta transformation.
 722   // In case that the TGeo was not initialized or not closed,
 723   // returns false and the object parameters are not set.
 724   //
 725   if (!gGeoManager || !gGeoManager->IsClosed()) {
 726     AliError("Can't set the local alignment object parameters! gGeoManager doesn't exist or it is still open!");
 727     return kFALSE;
 728   }
 729
 730   const char* symname = GetSymName();
 731   TGeoHMatrix gprime,gprimeinv;
 732   TGeoPhysicalNode* pn = 0;
 733   TGeoPNEntry* pne = gGeoManager->GetAlignableEntry(symname);
 734   if(pne)
 735   {
 736     pn = pne->GetPhysicalNode();
 737     if(pn){
 738       if (pn->IsAligned())
 739         AliWarning(Form("Volume %s has been misaligned already!",symname));
 740       gprime = *pn->GetMatrix();
 741     }else{
 742       gprime = pne->GetGlobalOrig();
 743     }
 744   }else{
 745     AliWarning(Form("The symbolic volume name %s does not correspond to a physical entry. Using it as volume path!",symname));
 746     if(!gGeoManager->cd(symname)) {
 747       AliError(Form("Volume name or path %s not valid!",symname));
 748       return kFALSE;
 749     }
 750     gprime = *gGeoManager->GetCurrentMatrix();
 751   }
 752
 753   TGeoHMatrix m1; // the TGeoHMatrix copy of the local delta "m"
 754   const Double_t *tr = m.GetTranslation();
 755   m1.SetTranslation(tr);
 756   const Double_t* rot = m.GetRotationMatrix();
 757   m1.SetRotation(rot);
 758
 759   gprimeinv = gprime.Inverse();
 760   m1.Multiply(&gprimeinv);
 761   m1.MultiplyLeft(&gprime);
 762
 763   return SetMatrix(m1);
 764 }
 765
 766 //_____________________________________________________________________________
 767 Bool_t AliAlignObj::SetMatrix(const TGeoMatrix& m)
 768 {
 769   // Set the global delta transformation by passing the TGeoMatrix
 770   // for it
 771   //
 772   SetTranslation(m);
 773   return SetRotation(m);
 774 }
 775
 776 //_____________________________________________________________________________
 777 Bool_t AliAlignObj::GetLocalPars(Double_t transl[], Double_t angles[]) const
 778 {
 779   // Get the translations and angles (in degrees) expressing the
 780   // local delta transformation.
 781   // In case that the TGeo was not initialized or not closed,
 782   // returns false and the object parameters are not set.
 783   //
 784   if(!GetLocalTranslation(transl)) return kFALSE;
 785   return GetLocalAngles(angles);
 786 }
 787
 788 //_____________________________________________________________________________
 789 Bool_t AliAlignObj::GetLocalTranslation(Double_t* tr) const
 790 {
 791   // Get the 3 shifts giving the translational part of the local
 792   // delta transformation.
 793   // In case that the TGeo was not initialized or not closed,
 794   // returns false and the object parameters are not set.
 795   //
 796   TGeoHMatrix ml;
 797   if(!GetLocalMatrix(ml)) return kFALSE;
 798   const Double_t* transl;
 799   transl = ml.GetTranslation();
 800   tr[0]=transl[0];
 801   tr[1]=transl[1];
 802   tr[2]=transl[2];
 803   return kTRUE;
 804 }
 805
 806 //_____________________________________________________________________________
 807 Bool_t AliAlignObj::GetLocalAngles(Double_t* angles) const
 808 {
 809   // Get the 3 angles giving the rotational part of the local
 810   // delta transformation.
 811   // In case that the TGeo was not initialized or not closed,
 812   // returns false and the object parameters are not set.
 813   //
 814   TGeoHMatrix ml;
 815   if(!GetLocalMatrix(ml)) return kFALSE;
 816   const Double_t *rot = ml.GetRotationMatrix();
 817   return MatrixToAngles(rot,angles);
 818 }
 819
 820 //_____________________________________________________________________________
 821 Bool_t AliAlignObj::GetLocalMatrix(TGeoHMatrix& m) const
 822 {
 823   // Get the matrix for the local delta transformation.
 824   // In case that the TGeo was not initialized or not closed,
 825   // returns false and the object parameters are not set.
 826   //
 827   if (!gGeoManager || !gGeoManager->IsClosed()) {
 828     AliError("Can't get the local alignment object parameters! gGeoManager doesn't exist or it is still open!");
 829     return kFALSE;
 830   }
 831
 832   const char* symname = GetSymName();
 833   TGeoPhysicalNode* node;
 834   TGeoPNEntry* pne = gGeoManager->GetAlignableEntry(symname);
 835   if(pne){
 836     if(!pne->GetPhysicalNode()){
 837       node = gGeoManager->MakeAlignablePN(pne);
 838     }else{
 839       node = pne->GetPhysicalNode();
 840     }
 841   }else{
 842     AliWarning(Form("The symbolic volume name %s does not correspond to a physical entry. Using it as volume path!",symname));
 843     node = (TGeoPhysicalNode*) gGeoManager->MakePhysicalNode(symname);
 844   }
 845
 846   if (!node) {
 847     AliError(Form("Volume name or path %s not valid!",symname));
 848     return kFALSE;
 849   }
 850 //  if (node->IsAligned())
 851 //    AliWarning(Form("Volume %s has been misaligned already!",symname));
 852
 853   GetMatrix(m);
 854   TGeoHMatrix gprime,gprimeinv;
 855   gprime = *node->GetMatrix();
 856   gprimeinv = gprime.Inverse();
 857   m.Multiply(&gprime);
 858   m.MultiplyLeft(&gprimeinv);
 859
 860   return kTRUE;
 861 }
 862
 863 //_____________________________________________________________________________
 864 Bool_t AliAlignObj::ApplyToGeometry(Bool_t ovlpcheck)
 865 {
 866   // Apply the current alignment object to the TGeo geometry
 867   // This method returns FALSE if the symname of the object was not
 868   // valid neither to get a TGeoPEntry nor as a volume path
 869   //
 870   if (!gGeoManager || !gGeoManager->IsClosed()) {
 871     AliError("Can't apply the alignment object! gGeoManager doesn't exist or it is still open!");
 872     return kFALSE;
 873   }
 874
 875   if (gGeoManager->IsLocked()){
 876     AliError("Can't apply the alignment object! Geometry is locked!");
 877     return kFALSE;
 878   }
 879
 880   const char* symname = GetSymName();
 881   const char* path;
 882   TGeoPhysicalNode* node;
 883   TGeoPNEntry* pne = gGeoManager->GetAlignableEntry(symname);
 884   if(pne){
 885     path = pne->GetTitle();
 886     node = gGeoManager->MakeAlignablePN(pne);
 887   }else{
 888     AliDebug(1,Form("The symbolic volume name %s does not correspond to a physical entry. Using it as a volume path!",symname));
 889     path=symname;
 890     if (!gGeoManager->CheckPath(path)) {
 891       AliDebug(1,Form("Volume path %s not valid!",path));
 892       return kFALSE;
 893     }
 894     if (gGeoManager->GetListOfPhysicalNodes()->FindObject(path)) {
 895       AliError(Form("Volume %s has been misaligned already!",path));
 896       return kFALSE;
 897     }
 898     node = (TGeoPhysicalNode*) gGeoManager->MakePhysicalNode(path);
 899   }
 900
 901   if (!node) {
 902     AliError(Form("Volume path %s not valid!",path));
 903     return kFALSE;
 904   }
 905
 906   //  Double_t threshold = 0.001;
 907
 908   TGeoHMatrix align,gprime;
 909   gprime = *node->GetMatrix();
 910   GetMatrix(align);
 911   gprime.MultiplyLeft(&align);
 912   TGeoHMatrix *ginv = new TGeoHMatrix;
 913   TGeoHMatrix *g = node->GetMatrix(node->GetLevel()-1);
 914   *ginv = g->Inverse();
 915   *ginv *= gprime;
 916   AliGeomManager::ELayerID layerId; // unique identity for layer in the alobj
 917   Int_t modId; // unique identity for volume inside layer in the alobj
 918   GetVolUID(layerId, modId);
 919   AliDebug(2,Form("Aligning volume %s of detector layer %d with local ID %d",symname,layerId,modId));
 920   if(ovlpcheck){
 921     node->Align(ginv,0,kTRUE); //(trunk of root takes threshold as additional argument)
 922   }else{
 923     node->Align(ginv,0,kFALSE);
 924   }
 925   if(ovlpcheck)
 926   {
 927     TObjArray* ovlpArray =  gGeoManager->GetListOfOverlaps();
 928     Int_t nOvlp = ovlpArray->GetEntriesFast();
 929     if(nOvlp)
 930     {
 931       AliInfo(Form("Misalignment of node %s generated the following overlaps/extrusions:",node->GetName()));
 932       for(Int_t i=0; i<nOvlp; i++)
 933         ((TGeoOverlap*)ovlpArray->UncheckedAt(i))->PrintInfo();
 934     }
 935   }
 936
 937   return kTRUE;
 938 }
 939
 940