HLT/TPCLib/AliHLTTPCClusterTransformation.cxx

   1 // $Id: AliHLTTPCClusterTransformation.cxx 41244 2010-05-14 08:13:35Z kkanaki $
   2
   3 //**************************************************************************
   4 //* This file is property of and copyright by the ALICE HLT Project        *
   5 //* ALICE Experiment at CERN, All rights reserved.                         *
   6 //*                                                                        *
   7 //* Primary Authors: Kalliopi Kanaki <Kalliopi.Kanaki@ift.uib.no>          *
   8 //*                  for The ALICE HLT Project.                            *
   9 //*                                                                        *
  10 //* Permission to use, copy, modify and distribute this software and its   *
  11 //* documentation strictly for non-commercial purposes is hereby granted   *
  12 //* without fee, provided that the above copyright notice appears in all   *
  13 //* copies and that both the copyright notice and this permission notice   *
  14 //* appear in the supporting documentation. The authors make no claims     *
  15 //* about the suitability of this software for any purpose. It is          *
  16 //* provided "as is" without express or implied warranty.                  *
  17 //**************************************************************************
  18
  19 /** @file   AliHLTTPCClusterTransformation.cxx
  20     @author Kalliopi Kanaki, Sergey Gorbubnov
  21     @date
  22     @brief
  23 */
  24
  25
  26 #include "AliHLTTPCClusterTransformation.h"
  27 #include "AliHLTTPCTransform.h"
  28
  29 #include "AliTPCcalibDB.h"
  30 #include "AliTPCTransform.h"
  31 #include "AliTPCParam.h"
  32 #include "AliTPCRecoParam.h"
  33 #include "AliGeomManager.h"
  34 #include <iostream>
  35 #include <iomanip>
  36
  37 using namespace std;
  38
  39 ClassImp(AliHLTTPCClusterTransformation) //ROOT macro for the implementation of ROOT specific class methods
  40
  41 AliHLTTPCClusterTransformation::AliHLTTPCClusterTransformation()
  42 :
  43   fOfflineTransform(NULL),
  44   fOfflineTPCParam( NULL ),
  45   fLastSector(-1)
  46 {
  47   // see header file for class documentation
  48   // or
  49   // refer to README to build package
  50   // or
  51   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  52   fAliT[0] = 0.;
  53   fAliT[1] = 0.;
  54   fAliT[2] = 0.;
  55   SetRotationMatrix();
  56 }
  57 AliHLTTPCClusterTransformation::~AliHLTTPCClusterTransformation()
  58 {
  59   // see header file for class documentation
  60   //delete fOfflineTransform;
  61 }
  62
  63
  64 int  AliHLTTPCClusterTransformation::Init( double FieldBz, UInt_t TimeStamp )
  65 {
  66   // Initialisation
  67
  68   //delete fOfflineTransform;
  69   fOfflineTransform = 0;
  70   fOfflineTPCParam = 0;
  71
  72   AliTPCcalibDB* pCalib=AliTPCcalibDB::Instance();
  73
  74   if(!pCalib ) return -1;
  75
  76   pCalib->SetExBField(FieldBz);
  77
  78   if(!AliGeomManager::GetGeometry()){
  79      AliGeomManager::LoadGeometry();
  80   }
  81
  82   if( !pCalib->GetTransform() ) return -2;
  83
  84   //fOfflineTransform = new AliTPCTransform (*pCalib->GetTransform());
  85   fOfflineTransform = pCalib->GetTransform();
  86   fOfflineTransform->SetCurrentRecoParam( AliTPCRecoParam::GetHLTParam() );
  87   fOfflineTransform->SetCurrentTimeStamp( TimeStamp );
  88   fOfflineTPCParam = pCalib->GetParameters();
  89   if( !fOfflineTPCParam ) return -3;
  90
  91   fOfflineTPCParam->Update();
  92   fOfflineTPCParam->ReadGeoMatrices();
  93
  94   fLastSector = -1;
  95
  96   fAliT[0] = 0.;
  97   fAliT[1] = 0.;
  98   fAliT[2] = 0.;
  99   SetRotationMatrix();
 100
 101   return 0;
 102 }
 103
 104
 105 void AliHLTTPCClusterTransformation::SetCurrentTimeStamp( UInt_t TimeStamp )
 106 {
 107   // Set the current time stamp
 108   if( fOfflineTransform ) fOfflineTransform->SetCurrentTimeStamp( TimeStamp );
 109 }
 110
 111
 112 int  AliHLTTPCClusterTransformation::Transform( int Slice, int Row, float Pad, float Time, float XYZ[] )
 113 {
 114   // Convert row, pad, time to X Y Z
 115
 116   Int_t sector=-99, thisrow=-99;
 117   AliHLTTPCTransform::Slice2Sector( Slice, Row, sector, thisrow);
 118
 119   if( !fOfflineTransform ){
 120     AliHLTTPCTransform::Raw2Local( XYZ, sector, thisrow, Pad, Time);
 121     if(Slice>17) XYZ[1]= - XYZ[1];
 122     return 0;
 123   }
 124
 125   Int_t iSector[1]= {sector};
 126   Double_t x[3] = { thisrow, Pad, Time };
 127   fOfflineTransform->Transform(x,iSector,0,1);
 128
 129   if( sector!= fLastSector ){
 130     if( fOfflineTPCParam && sector<fOfflineTPCParam->GetNSector() ){
 131       TGeoHMatrix  *alignment = fOfflineTPCParam->GetClusterMatrix( sector );
 132       if ( alignment ){
 133         const Double_t *tr = alignment->GetTranslation();
 134         const Double_t *rot = alignment->GetRotationMatrix();
 135         if( tr && rot ){
 136           for( int i=0; i<3; i++ ) fAliT[i] = tr[i];
 137           for( int i=0; i<9; i++ ) fAliR[i] = rot[i];
 138         }
 139       }
 140     } else {
 141       fAliT[0] = 0.;
 142       fAliT[1] = 0.;
 143       fAliT[2] = 0.;
 144       for( int i=0; i<9; i++ ) fAliR[i] = 0;
 145       fAliR[0] = 1.;
 146       fAliR[4] = 1.;
 147       fAliR[8] = 1.;
 148     }
 149     fLastSector = sector;
 150   }
 151   // alignment->LocalToMaster( x, y);
 152
 153   XYZ[0] = fAliT[0] + x[0]*fAliR[0] + x[1]*fAliR[1] + x[2]*fAliR[2];
 154   XYZ[1] = fAliT[1] + x[0]*fAliR[3] + x[1]*fAliR[4] + x[2]*fAliR[5];
 155   XYZ[2] = fAliT[2] + x[0]*fAliR[6] + x[1]*fAliR[7] + x[2]*fAliR[8];
 156
 157   return 0;
 158 }
 159
 160 int  AliHLTTPCClusterTransformation::ReverseAlignment( float XYZ[], int slice, int padrow)
 161 {
 162   // reverse the alignment correction
 163   Int_t sector=-99, thisrow=-99;
 164   AliHLTTPCTransform::Slice2Sector( slice, padrow, sector, thisrow);
 165   if( sector!= fLastSector ){
 166     if( fOfflineTPCParam && sector<fOfflineTPCParam->GetNSector() ){
 167       TGeoHMatrix  *alignment = fOfflineTPCParam->GetClusterMatrix( sector );
 168       if ( alignment ){
 169         const Double_t *tr = alignment->GetTranslation();
 170         const Double_t *rot = alignment->GetRotationMatrix();
 171         if(tr){
 172           for( int i=0; i<3; i++ ) fAliT[i] = tr[i];
 173         }
 174         SetRotationMatrix(rot, true);
 175       }
 176     } else {
 177       fAliT[0] = 0.;
 178       fAliT[1] = 0.;
 179       fAliT[2] = 0.;
 180       SetRotationMatrix(NULL, true);
 181     }
 182     fLastSector = sector;
 183   }
 184
 185   // correct for alignment: translation
 186   float xyz[3];
 187   xyz[0] = XYZ[0] - fAliT[0];
 188   xyz[1] = XYZ[1] - fAliT[1];
 189   xyz[2] = XYZ[2] - fAliT[2];
 190
 191   // correct for alignment: rotation
 192   XYZ[0]=xyz[0]*fAdjR[0] + xyz[1]*fAdjR[1] + xyz[2]*fAdjR[2];
 193   XYZ[1]=xyz[0]*fAdjR[3] + xyz[1]*fAdjR[4] + xyz[2]*fAdjR[5];
 194   XYZ[2]=xyz[0]*fAdjR[6] + xyz[1]*fAdjR[7] + xyz[2]*fAdjR[8];
 195
 196   return 0;
 197 }
 198
 199 void AliHLTTPCClusterTransformation::SetRotationMatrix(const Double_t *rot, bool bCalcAdjugate)
 200 {
 201   // set the rotation matrix and calculate the adjugate if requested
 202   if (rot) {
 203     for( int i=0; i<9; i++ ) fAliR[i] = rot[i];
 204     if (bCalcAdjugate) {
 205       CalcAdjugateRotation();
 206     }
 207     return;
 208   }
 209   for( int i=0; i<9; i++ ) {fAliR[i] = 0; fAdjR[i] = 0;}
 210   fAliR[0] = 1.;
 211   fAliR[4] = 1.;
 212   fAliR[8] = 1.;
 213   fAdjR[0] = 1.;
 214   fAdjR[4] = 1.;
 215   fAdjR[8] = 1.;
 216 }
 217
 218 bool AliHLTTPCClusterTransformation::CalcAdjugateRotation(bool bCheck)
 219 {
 220   // check rotation matrix and adjugate for consistency
 221   fAdjR[0]= fAliR[4]*fAliR[8]-fAliR[5]*fAliR[7];
 222   fAdjR[1]= fAliR[5]*fAliR[6]-fAliR[3]*fAliR[8];
 223   fAdjR[2]= fAliR[3]*fAliR[7]-fAliR[4]*fAliR[6];
 224
 225   fAdjR[3]= fAliR[2]*fAliR[7]-fAliR[1]*fAliR[8];
 226   fAdjR[4]= fAliR[0]*fAliR[8]-fAliR[2]*fAliR[6];
 227   fAdjR[5]= fAliR[2]*fAliR[6]-fAliR[0]*fAliR[7];
 228
 229   fAdjR[6]= fAliR[1]*fAliR[5]-fAliR[2]*fAliR[4];
 230   fAdjR[7]= fAliR[2]*fAliR[3]-fAliR[0]*fAliR[5];
 231   fAdjR[8]= fAliR[0]*fAliR[4]-fAliR[1]*fAliR[3];
 232
 233   if (bCheck) {
 234     for (int r=0; r<3; r++) {
 235       for (int c=0; c<3; c++) {
 236         float a=0.;
 237         float expected=0.;
 238         if (r==c) expected=1.;
 239         for (int i=0; i<3; i++) {
 240           a+=fAliR[3*r+i]*fAdjR[c+(3*i)];
 241         }
 242         if (TMath::Abs(a-expected)>0.00001) {
 243           std::cout << "inconsistent adjugate at " << r << c << ": " << a << " " << expected << std::endl;
 244           return false;
 245         }
 246       }
 247     }
 248   }
 249   return true;
 250 }
 251
 252 void AliHLTTPCClusterTransformation::Print(const char* /*option*/) const
 253 {
 254   // print info
 255   ios::fmtflags coutflags=std::cout.flags(); // backup cout status flags
 256   std::cout << "AliHLTTPCClusterTransformation for sector " << fLastSector << std::endl;
 257
 258   std::cout.setf(ios_base::showpos|ios_base::showpos|ios::right);
 259   std::cout << "  translation: " << std::endl;
 260   int r=0;
 261   for (r=0; r<3; r++) {
 262     std::cout << setw(7) << fixed << setprecision(2);
 263     cout << "  " << fAliT[r] << std::endl;
 264   }
 265   std::cout << "  rotation and adjugated rotation: " << std::endl;
 266   for (r=0; r<3; r++) {
 267     int c=0;
 268     std::cout << setw(7) << fixed << setprecision(2);
 269     for (c=0; c<3; c++) std::cout << "  " << fAliR[3*r+c];
 270     std::cout << "      ";
 271     for (c=0; c<3; c++) std::cout << "  " << fAdjR[3*r+c];
 272     std::cout << endl;
 273   }
 274   std::cout.flags(coutflags); // restore the original flags
 275 }