HLT/TPCLib/transform/AliHLTTPCFastTransform.cxx

   1 //**************************************************************************
   2 //* This file is property of and copyright by the ALICE HLT Project        *
   3 //* ALICE Experiment at CERN, All rights reserved.                         *
   4 //*                                                                        *
   5 //* Primary Authors: Sergey Gorbunov <sergey.gorbunov@cern.ch>             *
   6 //*                  for The ALICE HLT Project.                            *
   7 //*                                                                        *
   8 //* Permission to use, copy, modify and distribute this software and its   *
   9 //* documentation strictly for non-commercial purposes is hereby granted   *
  10 //* without fee, provided that the above copyright notice appears in all   *
  11 //* copies and that both the copyright notice and this permission notice   *
  12 //* appear in the supporting documentation. The authors make no claims     *
  13 //* about the suitability of this software for any purpose. It is          *
  14 //* provided "as is" without express or implied warranty.                  *
  15 //**************************************************************************
  16
  17 /** @file   AliHLTTPCFastTransform.cxx
  18     @author Sergey Gorbubnov
  19     @date
  20     @brief
  21 */
  22
  23
  24 #include "AliHLTTPCFastTransform.h"
  25 #include "AliTPCTransform.h"
  26 #include "AliTPCParam.h"
  27 #include "AliTPCcalibDB.h"
  28
  29 #include <iostream>
  30 #include <iomanip>
  31
  32 using namespace std;
  33
  34 ClassImp(AliHLTTPCFastTransform); //ROOT macro for the implementation of ROOT specific class methods
  35
  36
  37 AliHLTTPCFastTransform::AliHLTTPCFastTransform()
  38 :
  39   fError(),
  40   fOrigTransform(0),
  41   fLastTimeStamp(-1),
  42   fLastTimeBin(600),
  43   fTimeBorder1(100),
  44   fTimeBorder2(500),
  45   fAlignment(NULL)
  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   for( Int_t i=0; i<fkNSec; i++)
  53     for( Int_t j=0; j<fkNRows; j++ ) fRows[i][j] = NULL;
  54 }
  55
  56 AliHLTTPCFastTransform::~AliHLTTPCFastTransform()
  57 {
  58   // see header file for class documentation
  59   DeInit();
  60 }
  61
  62 void  AliHLTTPCFastTransform::DeInit()
  63 {
  64   // Deinitialisation
  65
  66   for( Int_t i=0; i<fkNSec; i++){
  67     for( Int_t j=0; j<fkNRows; j++ ){
  68       delete fRows[i][j];
  69       fRows[i][j] = 0;
  70     }
  71   }
  72   fOrigTransform = NULL;
  73   fLastTimeStamp = -1;
  74   delete[] fAlignment;
  75   fAlignment = NULL;
  76   fError = "";
  77 }
  78
  79
  80 Int_t  AliHLTTPCFastTransform::Init( AliTPCTransform *transform, Long_t TimeStamp )
  81 {
  82   // Initialisation
  83
  84   DeInit();
  85
  86   AliTPCcalibDB* pCalib=AliTPCcalibDB::Instance();
  87   if(!pCalib ) return Error( -1, "AliHLTTPCFastTransform::Init: No TPC calibration instance found");
  88
  89   AliTPCParam *tpcParam = pCalib->GetParameters();
  90   if( !tpcParam ) return Error( -2, "AliHLTTPCFastTransform::Init: No TPCParam object found");
  91
  92   if( !transform ) transform = pCalib->GetTransform();
  93   if( !transform ) return Error( -3, "AliHLTTPCFastTransform::Init: No TPC transformation found");
  94
  95   tpcParam->Update();
  96   tpcParam->ReadGeoMatrices();
  97
  98
  99   // at the moment initialise all the rows
 100
 101   int nSec = tpcParam->GetNSector();
 102   if( nSec>fkNSec ) nSec = fkNSec;
 103
 104   for( Int_t i=0; i<nSec; i++ ){
 105     int nRows = tpcParam->GetNRow(i);
 106     if( nRows>fkNRows ) nRows = fkNRows;
 107     for( int j=0; j<nRows; j++){
 108       if( !fRows[i][j] ) fRows[i][j] = new AliHLTTPCFastTransform::AliRowTransform;
 109       if( !fRows[i][j] ) return Error( -4, "AliHLTTPCFastTransform::Init: Not enough memory");
 110     }
 111   }
 112
 113   const AliTPCRecoParam *rec = transform->GetCurrentRecoParam();
 114
 115   if( !rec ) return Error( -5, "AliHLTTPCFastTransform::Init: No TPC Reco Param set in transformation");
 116
 117   fOrigTransform = transform;
 118
 119   if( rec->GetUseSectorAlignment() && (!pCalib->HasAlignmentOCDB()) ){
 120
 121     fAlignment = new Float_t [fkNSec*21];
 122     for( Int_t iSec=0; iSec<fkNSec; iSec++ ){
 123       Float_t *t = fAlignment + iSec*21, *r = t+3, *v = t+12;
 124       for( int i=0; i<21; i++ ) t[i]=0.f;
 125       r[0] = r[4] = r[8] = 1.f;
 126       v[0] = v[4] = v[8] = 1.f;
 127     }
 128
 129     for( Int_t iSec=0; iSec<nSec; iSec++ ){
 130       Float_t *t = fAlignment + iSec*21, *r = t+3, *v = t+12;
 131       TGeoHMatrix  *alignment = tpcParam->GetClusterMatrix( iSec );
 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++ ) t[i] = tr[i];
 137           for( int i=0; i<9; i++ ) r[i] = rot[i];
 138           CalcAdjugateRotation(r,v,1);
 139         }
 140       }
 141     }
 142   }
 143
 144   return SetCurrentTimeStamp( TimeStamp );
 145 }
 146
 147
 148 bool AliHLTTPCFastTransform::CalcAdjugateRotation(const Float_t *mA, Float_t *mB, bool bCheck)
 149 {
 150   // check rotation matrix and adjugate for consistency
 151
 152   mB[0]= mA[4]*mA[8]-mA[5]*mA[7];
 153   mB[1]= mA[5]*mA[6]-mA[3]*mA[8];
 154   mB[2]= mA[3]*mA[7]-mA[4]*mA[6];
 155
 156   mB[3]= mA[2]*mA[7]-mA[1]*mA[8];
 157   mB[4]= mA[0]*mA[8]-mA[2]*mA[6];
 158   mB[5]= mA[2]*mA[6]-mA[0]*mA[7];
 159
 160   mB[6]= mA[1]*mA[5]-mA[2]*mA[4];
 161   mB[7]= mA[2]*mA[3]-mA[0]*mA[5];
 162   mB[8]= mA[0]*mA[4]-mA[1]*mA[3];
 163
 164   if (bCheck) {
 165     for (int r=0; r<3; r++) {
 166       for (int c=0; c<3; c++) {
 167         float a=0.;
 168         float expected=0.;
 169         if (r==c) expected=1.;
 170         for (int i=0; i<3; i++) {
 171           a+=mA[3*r+i]*mB[c+(3*i)];
 172         }
 173         if (TMath::Abs(a-expected)>0.00001) {
 174           std::cout << "inconsistent adjugate at " << r << c << ": " << a << " " << expected << std::endl;
 175           for( int i=0; i<9; i++ ) mB[i] = 0;
 176           mB[0] = mB[4] = mB[8] = 1;
 177           return false;
 178         }
 179       }
 180     }
 181   }
 182   return true;
 183 }
 184
 185
 186
 187 Int_t AliHLTTPCFastTransform::SetCurrentTimeStamp( Long_t TimeStamp )
 188 {
 189   // Set the current time stamp
 190
 191
 192   Long_t lastTS = fLastTimeStamp;
 193   fLastTimeStamp = -1; // deinitialise
 194
 195   if( !fOrigTransform ) return Error( -1, "AliHLTTPCFastTransform::SetCurrentTimeStamp: TPC transformation has not been set properly");
 196
 197   AliTPCcalibDB* pCalib=AliTPCcalibDB::Instance();
 198   if(!pCalib ) return Error( -2, "AliHLTTPCFastTransform::SetCurrentTimeStamp: No TPC calibration found");
 199
 200   AliTPCParam *tpcParam = pCalib->GetParameters();
 201   if( !tpcParam ) return  Error( -3, "AliHLTTPCFastTransform::SetCurrentTimeStamp: No TPCParam object found");
 202
 203
 204   if( TimeStamp<0  ) return 0;
 205
 206   fLastTimeStamp = lastTS;
 207
 208   if( fLastTimeStamp>=0 && TMath::Abs(fLastTimeStamp - TimeStamp ) <60 ) return 0;
 209
 210
 211   fOrigTransform->SetCurrentTimeStamp( static_cast<UInt_t>(TimeStamp) );
 212   fLastTimeStamp = TimeStamp;
 213
 214   // find last calibrated time bin
 215
 216   Int_t nTimeBins = tpcParam->GetMaxTBin();
 217   Int_t is[]={0};
 218   bool sign = 0;
 219   for( fLastTimeBin=0; fLastTimeBin<nTimeBins; fLastTimeBin++){
 220     Double_t xx[]={0,0,fLastTimeBin};
 221     fOrigTransform->Transform(xx,is,0,1);
 222     bool s = (xx[2]>=0);
 223     if( fLastTimeBin==0 ) sign = s;
 224     else if( sign!=s ){
 225       fLastTimeBin--;
 226       break;
 227     }
 228   }
 229   fTimeBorder1 = 60;
 230   fTimeBorder2 = fLastTimeBin - 100;
 231
 232   int nSec = tpcParam->GetNSector();
 233   if( nSec>fkNSec ) nSec = fkNSec;
 234
 235   for( Int_t i=0; i<nSec; i++ ){
 236     int nRows = tpcParam->GetNRow(i);
 237     if( nRows>fkNRows ) nRows = fkNRows;
 238     for( int j=0; j<nRows; j++){
 239       if( fRows[i][j] ){
 240         int err = InitRow(i,j);
 241         if( err!=0 ) return err;
 242       }
 243     }
 244   }
 245   return 0;
 246 }
 247
 248
 249 Int_t AliHLTTPCFastTransform::InitRow( Int_t iSector, Int_t iRow )
 250 {
 251   // see header file for class documentation
 252
 253   AliTPCcalibDB* pCalib=AliTPCcalibDB::Instance();
 254
 255   if( iSector<0 || iSector>=fkNSec || iRow<0 || iRow>=fkNRows || !fOrigTransform || (fLastTimeStamp<0) ||
 256       !fRows[iSector][iRow] || !pCalib || !pCalib->GetParameters() ){
 257     return Error( -1, "AliHLTTPCFastTransform::InitRow: Internal error");
 258   }
 259
 260   AliTPCParam *tpcParam = pCalib->GetParameters();
 261
 262   Int_t nPads = tpcParam->GetNPads(iSector,iRow);
 263
 264   if( nPads<2 ){
 265     return Error( -2, Form("AliHLTTPCFastTransform::InitRow: Wrong NPads=%d for sector %d row %d",nPads,iSector,iRow ));
 266   }
 267
 268   fRows[iSector][iRow]->fSpline[0].Init(         0.5, nPads-1+0.5, 15, 0,  fTimeBorder1,         5);
 269   fRows[iSector][iRow]->fSpline[1].Init(         0.5, nPads-1+0.5, 15, fTimeBorder1, fTimeBorder2,         10);
 270   fRows[iSector][iRow]->fSpline[2].Init(         0.5, nPads-1+0.5, 15, fTimeBorder2, fLastTimeBin,         5);
 271
 272   for( Int_t i=0; i<3; i++){
 273     Int_t is[]={iSector};
 274     for( Int_t j=0; j<fRows[iSector][iRow]->fSpline[i].GetNPoints(); j++){
 275       Float_t pad, time;
 276       fRows[iSector][iRow]->fSpline[i].GetAB(j,pad,time);
 277       Double_t xx[]={iRow,pad,time};
 278       fOrigTransform->Transform(xx,is,0,1);
 279       fRows[iSector][iRow]->fSpline[i].Fill(j,xx);
 280     }
 281     fRows[iSector][iRow]->fSpline[i].Consolidate();
 282   }
 283   return 0;
 284 }
 285
 286
 287
 288 Int_t AliHLTTPCFastTransform::GetRowSize( Int_t iSec, Int_t iRow ) const
 289 {
 290   // see header file for class documentation
 291   Int_t s = sizeof(AliHLTTPCFastTransform::AliRowTransform);
 292   if( fRows[iSec][iRow] ) for( Int_t i=0; i<3; i++) s+=fRows[iSec][iRow]->fSpline[i].GetMapSize();
 293   return s;
 294 }
 295
 296 Int_t AliHLTTPCFastTransform::GetSize() const
 297 {
 298   // see header file for class documentation
 299   Int_t s = sizeof(AliHLTTPCFastTransform);
 300   for( Int_t i=0; i<fkNSec; i++ )
 301     for( Int_t j=0; j<fkNRows; j++ ) if( fRows[i][j] ){
 302         s+= sizeof(AliHLTTPCFastTransform::AliRowTransform);
 303         for( Int_t k=0; k<3; k++) fRows[i][j]->fSpline[k].GetMapSize();
 304       }
 305   return s;
 306 }
 307
 308
 309 void AliHLTTPCFastTransform::Print(const char* /*option*/) const
 310 {
 311   // print info
 312   ios::fmtflags coutflags=std::cout.flags(); // backup cout status flags
 313
 314   if( !fAlignment ){
 315     std::cout << "AliHLTTPCFastTransform: no alignment transformation";
 316   } else {
 317     for( int iSec=0; iSec<fkNSec; iSec++ ){
 318       std::cout << "AliHLTTPCClusterTransformation for sector " << iSec << std::endl;
 319
 320       const Float_t *mT = fAlignment + iSec*21;
 321       const Float_t *mR = mT + 3;
 322       const Float_t *mV = mT + 3+9;
 323
 324       std::cout.setf(ios_base::showpos|ios_base::showpos|ios::right);
 325       std::cout << "  translation: " << std::endl;
 326       int r=0;
 327       for (r=0; r<3; r++) {
 328         std::cout << setw(7) << fixed << setprecision(2);
 329         cout << "  " << mT[r] << std::endl;
 330       }
 331       std::cout << "  rotation and adjugated rotation: " << std::endl;
 332       for (r=0; r<3; r++) {
 333         int c=0;
 334         std::cout << setw(7) << fixed << setprecision(2);
 335         for (c=0; c<3; c++) std::cout << "  " << mR[3*r+c];
 336         std::cout << "      ";
 337         for (c=0; c<3; c++) std::cout << "  " << mV[3*r+c];
 338         std::cout << endl;
 339       }
 340     }
 341   }
 342   std::cout.flags(coutflags); // restore the original flags
 343 }