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