TPC/Base/AliTPCCorrectionLookupTable.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 #include <TMath.h>
  17 #include <TMatrixF.h>
  18 #include <TStopwatch.h>
  19 #include <TString.h>
  20 #include <TFile.h>
  21 #include <TObjArray.h>
  22 #include <TSystem.h>
  23 #include <THashList.h>
  24 #include <TVector2.h>
  25 #include <TLinearFitter.h>
  26
  27 #include <AliLog.h>
  28 #include <AliTPCROC.h>
  29
  30 #include "AliTPCCorrection.h"
  31
  32 #include "AliTPCCorrectionLookupTable.h"
  33
  34 ClassImp(AliTPCCorrectionLookupTable)
  35
  36 //_________________________________________________________________________________________
  37 AliTPCCorrectionLookupTable::AliTPCCorrectionLookupTable()
  38 : AliTPCCorrection()
  39 , fNR(0)
  40 , fNPhi(0)
  41 , fNZ(0)
  42 , fCorrScaleFactor(-1)
  43 , fFillCorrection(kTRUE)
  44 , fLimitsR()
  45 , fLimitsPhi()
  46 , fLimitsZ()
  47 , fLookUpDxDist(0x0)
  48 , fLookUpDyDist(0x0)
  49 , fLookUpDzDist(0x0)
  50 , fLookUpDxCorr(0x0)
  51 , fLookUpDyCorr(0x0)
  52 , fLookUpDzCorr(0x0)
  53 {
  54   //
  55   //
  56   //
  57 }
  58 //_________________________________________________________________________________________
  59 AliTPCCorrectionLookupTable::~AliTPCCorrectionLookupTable()
  60 {
  61   //
  62   // dtor
  63   //
  64
  65   ResetTables();
  66 }
  67
  68 //_________________________________________________________________________________________
  69 void AliTPCCorrectionLookupTable::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
  70   //
  71   // Get interpolated Distortion
  72   //
  73
  74   GetInterpolation(x,roc,dx,fLookUpDxDist,fLookUpDyDist,fLookUpDzDist);
  75 }
  76
  77 //_________________________________________________________________________________________
  78 void AliTPCCorrectionLookupTable::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
  79   //
  80   // Get interplolated correction
  81   //
  82   GetInterpolation(x,roc,dx,fLookUpDxCorr,fLookUpDyCorr,fLookUpDzCorr);
  83
  84   if (fCorrScaleFactor>0) {
  85     dx[0]*=fCorrScaleFactor;
  86     dx[1]*=fCorrScaleFactor;
  87     dx[2]*=fCorrScaleFactor;
  88   }
  89 }
  90
  91 //_________________________________________________________________________________________
  92 void AliTPCCorrectionLookupTable::GetInterpolation(const Float_t x[],const Short_t roc,Float_t dx[],
  93                                                    TMatrixF **mDx, TMatrixF **mDy, TMatrixF **mDz)
  94 {
  95   //
  96   // Calculates the correction/distotring from a lookup table
  97   // type: 0 = correction
  98   //       1 = distortion
  99   //
 100
 101 //   Float_t typeSign=-1;
 102 //   if (type==1) typeSign=1;
 103
 104   Int_t   order     = 1 ;    // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
 105
 106   Double_t r, phi, z ;
 107   Int_t    sign;
 108
 109   r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
 110   phi    =  TMath::ATan2(x[1],x[0]) ;
 111   if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
 112   z      =  x[2] ;                                         // Create temporary copy of x[2]
 113
 114   if ( (roc%36) < 18 ) {
 115     sign =  1;       // (TPC A side)
 116   } else {
 117     sign = -1;       // (TPC C side)
 118   }
 119
 120   if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
 121   if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
 122
 123
 124   if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
 125     AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
 126
 127   // Get the Er and Ephi field integrals plus the integral over Z
 128     dx[0] = Interpolate3DTable(order, r, z, phi,
 129                                fNR, fNZ, fNPhi,
 130                                fLimitsR.GetMatrixArray(),
 131                                fLimitsZ.GetMatrixArray(),
 132                                fLimitsPhi.GetMatrixArray(),
 133                                mDx  );
 134     dx[1] = Interpolate3DTable(order, r, z, phi,
 135                                fNR, fNZ, fNPhi,
 136                                fLimitsR.GetMatrixArray(),
 137                                fLimitsZ.GetMatrixArray(),
 138                                fLimitsPhi.GetMatrixArray(),
 139                                mDy);
 140     dx[2] = Interpolate3DTable(order, r, z, phi,
 141                                fNR, fNZ, fNPhi,
 142                                fLimitsR.GetMatrixArray(),
 143                                fLimitsZ.GetMatrixArray(),
 144                                fLimitsPhi.GetMatrixArray(),
 145                                mDz   );
 146 }
 147
 148 //_________________________________________________________________________________________
 149 void AliTPCCorrectionLookupTable::CreateLookupTable(AliTPCCorrection &tpcCorr, Float_t stepSize/*=5.*/)
 150 {
 151   //
 152   // create lookup table for all phi,r,z bins
 153   //
 154
 155   if (fNR==0) {
 156     AliError("Limits are not set yet. Please use one of the Set..Limits functions first");
 157     return;
 158   }
 159
 160   TStopwatch s;
 161
 162   ResetTables();
 163   InitTables();
 164
 165   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 166     CreateLookupTablePhiBin(tpcCorr,iPhi,stepSize);
 167   }
 168
 169   s.Stop();
 170   AliInfo(Form("Required time for lookup table creation: %.2f (%.2f) sec. real (cpu)",s.RealTime(),s.CpuTime()));
 171 }
 172
 173 //_________________________________________________________________________________________
 174 void AliTPCCorrectionLookupTable::CreateLookupTableSinglePhi(AliTPCCorrection &tpcCorr, Int_t iPhi, Float_t stepSize)
 175 {
 176   //
 177   // Lookup table for only one phi bin. Can be used for parallel processing
 178   //
 179
 180   if (fNR==0) {
 181     AliError("Limits are not set yet. Please use one of the Set..Limits functions first");
 182     return;
 183   }
 184
 185   TStopwatch s;
 186
 187   ResetTables();
 188   InitTableArrays();
 189   InitTablesPhiBin(iPhi);
 190   CreateLookupTablePhiBin(tpcCorr,iPhi,stepSize);
 191
 192   s.Stop();
 193   AliInfo(Form("Required time for lookup table creation: %.2f (%.2f) sec. real (cpu)",s.RealTime(),s.CpuTime()));
 194 }
 195
 196 //_________________________________________________________________________________________
 197 void AliTPCCorrectionLookupTable::CreateLookupTablePhiBin(AliTPCCorrection &tpcCorr, Int_t iPhi, Float_t stepSize)
 198 {
 199   //
 200   //
 201   //
 202
 203   if (iPhi<0||iPhi>=fNPhi) return;
 204
 205   const Float_t delta=stepSize; // 5cm
 206   Float_t x[3]={0.,0.,0.};
 207   Float_t dx[3]={0.,0.,0.};
 208
 209   Double_t phi=fLimitsPhi(iPhi);
 210   //
 211   TMatrixF &mDxDist   = *fLookUpDxDist[iPhi];
 212   TMatrixF &mDyDist   = *fLookUpDyDist[iPhi];
 213   TMatrixF &mDzDist   = *fLookUpDzDist[iPhi];
 214   //
 215   TMatrixF &mDxCorr   = *fLookUpDxCorr[iPhi];
 216   TMatrixF &mDyCorr   = *fLookUpDyCorr[iPhi];
 217   TMatrixF &mDzCorr   = *fLookUpDzCorr[iPhi];
 218
 219   for (Int_t ir=0; ir<fNR; ++ir){
 220     Double_t r=fLimitsR(ir);
 221     x[0]=r * TMath::Cos(phi);
 222     x[1]=r * TMath::Sin(phi);
 223
 224     for (Int_t iz=0; iz<fNZ; ++iz){
 225       Double_t z=fLimitsZ(iz);
 226       x[2]=z;
 227       //TODO: change hardcoded value for r>133.?
 228       Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
 229       if (r>133.) roc+=36;
 230       if (z<0)    roc+=18;
 231
 232       if (delta>0)
 233         tpcCorr.GetDistortionIntegralDz(x,roc,dx,delta);
 234       else
 235         tpcCorr.GetDistortion(x,roc,dx);
 236       mDxDist(ir,iz)=dx[0];
 237       mDyDist(ir,iz)=dx[1];
 238       mDzDist(ir,iz)=dx[2];
 239
 240       if (fFillCorrection) {
 241         if (delta>0)
 242           tpcCorr.GetCorrectionIntegralDz(x,roc,dx,delta);
 243         else
 244           tpcCorr.GetCorrection(x,roc,dx);
 245         mDxCorr(ir,iz)=dx[0];
 246         mDyCorr(ir,iz)=dx[1];
 247         mDzCorr(ir,iz)=dx[2];
 248       }
 249     }
 250   }
 251
 252 }
 253
 254 //_________________________________________________________________________________________
 255 void AliTPCCorrectionLookupTable::InitTables()
 256 {
 257   //
 258   // Init all tables
 259   //
 260
 261   InitTableArrays();
 262   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 263     InitTablesPhiBin(iPhi);
 264   }
 265 }
 266
 267 //_________________________________________________________________________________________
 268 void AliTPCCorrectionLookupTable::CreateLookupTableFromResidualDistortion(THn &resDist)
 269 {
 270   //
 271   // create lookup table from residual distortions stored in a 3d histogram
 272   // assume dimensions are r, phi, z
 273   //
 274   if (fNR==0) {
 275     AliError("Limits are not set yet. Please use one of the Set..Limits functions first");
 276     return;
 277   }
 278
 279   ResetTables();
 280   InitTables();
 281
 282   Double_t x[3]={0.,0.,0.};
 283
 284   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 285     const Double_t phi=fLimitsPhi(iPhi);
 286     x[1]=phi;
 287     //
 288     TMatrixF &mDxDist   = *fLookUpDxDist[iPhi];
 289     TMatrixF &mDyDist   = *fLookUpDyDist[iPhi];
 290     TMatrixF &mDzDist   = *fLookUpDzDist[iPhi];
 291     //
 292     TMatrixF &mDxCorr   = *fLookUpDxCorr[iPhi];
 293     TMatrixF &mDyCorr   = *fLookUpDyCorr[iPhi];
 294     TMatrixF &mDzCorr   = *fLookUpDzCorr[iPhi];
 295
 296     for (Int_t ir=0; ir<fNR; ++ir){
 297       const Double_t r=fLimitsR(ir);
 298       x[0]=r;
 299
 300       for (Int_t iz=0; iz<fNZ; ++iz){
 301         const Double_t z=fLimitsZ(iz);
 302         x[2]=z;
 303
 304         const Double_t drphi = resDist.GetBinContent(resDist.GetBin(x));
 305         Double_t dx[3]={0.,drphi,0.};
 306
 307         // transform rphi distortions (local y, so dy') to a global distortion
 308         // assume no radial distortion (dx' = 0)
 309         // assume no residual distortion in z for the moment
 310         Double_t cs=TMath::Cos(phi), sn=TMath::Sin(phi), lx=dx[0];
 311         dx[0]=lx*cs - dx[1]*sn; dx[1]=lx*sn + dx[1]*cs;
 312
 313         mDxDist(ir,iz)=dx[0];
 314         mDyDist(ir,iz)=dx[1];
 315         mDzDist(ir,iz)=dx[2];
 316
 317         mDxCorr(ir,iz)=-dx[0];
 318         mDyCorr(ir,iz)=-dx[1];
 319         mDzCorr(ir,iz)=-dx[2];
 320       }
 321     }
 322   }
 323 }
 324
 325 //_________________________________________________________________________________________
 326 void AliTPCCorrectionLookupTable::CreateResidual(AliTPCCorrection *distortion, AliTPCCorrection* correction)
 327 {
 328   //
 329   // create lookup table from residual distortions calculated from distorted - correction
 330   //
 331
 332   ResetTables();
 333   InitTables();
 334
 335   Float_t x[3]={0.,0.,0.};
 336
 337   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 338     const Double_t phi=fLimitsPhi(iPhi);
 339     //
 340     TMatrixF &mDxDist   = *fLookUpDxDist[iPhi];
 341     TMatrixF &mDyDist   = *fLookUpDyDist[iPhi];
 342     TMatrixF &mDzDist   = *fLookUpDzDist[iPhi];
 343     //
 344     TMatrixF &mDxCorr   = *fLookUpDxCorr[iPhi];
 345     TMatrixF &mDyCorr   = *fLookUpDyCorr[iPhi];
 346     TMatrixF &mDzCorr   = *fLookUpDzCorr[iPhi];
 347
 348     for (Int_t ir=0; ir<fNR; ++ir){
 349       const Double_t r=fLimitsR(ir);
 350       x[0]=r * TMath::Cos(phi);
 351       x[1]=r * TMath::Sin(phi);
 352
 353       for (Int_t iz=0; iz<fNZ; ++iz){
 354         const Double_t z=fLimitsZ(iz);
 355         x[2]=z;
 356
 357         //original point
 358         Float_t xdc[3]={x[0], x[1], x[2]};
 359
 360         Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
 361         if (r>133.) roc+=36;
 362         if (z<0)    roc+=18;
 363
 364         //get residual distortion
 365         distortion->DistortPoint(xdc, roc);
 366         correction->CorrectPoint(xdc, roc);
 367         Float_t dx[3]={x[0]-xdc[0], x[1]-xdc[1], x[2]-xdc[2]};
 368
 369         mDxDist(ir,iz)=dx[0];
 370         mDyDist(ir,iz)=dx[1];
 371         mDzDist(ir,iz)=dx[2];
 372
 373         mDxCorr(ir,iz)=-dx[0];
 374         mDyCorr(ir,iz)=-dx[1];
 375         mDzCorr(ir,iz)=-dx[2];
 376       }
 377     }
 378   }
 379 }
 380
 381 //_________________________________________________________________________________________
 382 void AliTPCCorrectionLookupTable::InitTablesPhiBin(Int_t iPhi)
 383 {
 384   //
 385   //
 386   //
 387
 388   // check if already initialised
 389   if (iPhi<0||iPhi>=fNPhi) return;
 390   if (fLookUpDxCorr[iPhi]) return;
 391
 392   fLookUpDxCorr[iPhi] = new TMatrixF(fNR,fNZ);
 393   fLookUpDyCorr[iPhi] = new TMatrixF(fNR,fNZ);
 394   fLookUpDzCorr[iPhi] = new TMatrixF(fNR,fNZ);
 395
 396   fLookUpDxDist[iPhi] = new TMatrixF(fNR,fNZ);
 397   fLookUpDyDist[iPhi] = new TMatrixF(fNR,fNZ);
 398   fLookUpDzDist[iPhi] = new TMatrixF(fNR,fNZ);
 399
 400 }
 401 //_________________________________________________________________________________________
 402 void AliTPCCorrectionLookupTable::InitTableArrays()
 403 {
 404   //
 405   //
 406   //
 407
 408   // needs to be before the table creation to set the limits
 409   SetupDefaultLimits();
 410
 411   fLookUpDxCorr = new TMatrixF*[fNPhi];
 412   fLookUpDyCorr = new TMatrixF*[fNPhi];
 413   fLookUpDzCorr = new TMatrixF*[fNPhi];
 414
 415   fLookUpDxDist = new TMatrixF*[fNPhi];
 416   fLookUpDyDist = new TMatrixF*[fNPhi];
 417   fLookUpDzDist = new TMatrixF*[fNPhi];
 418
 419   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 420     fLookUpDxCorr[iPhi] = 0x0;
 421     fLookUpDyCorr[iPhi] = 0x0;
 422     fLookUpDzCorr[iPhi] = 0x0;
 423
 424     fLookUpDxDist[iPhi] = 0x0;
 425     fLookUpDyDist[iPhi] = 0x0;
 426     fLookUpDzDist[iPhi] = 0x0;
 427   }
 428 }
 429
 430 //_________________________________________________________________________________________
 431 void AliTPCCorrectionLookupTable::ResetTables()
 432 {
 433   //
 434   // Reset the lookup tables
 435   //
 436
 437   if (!fLookUpDxCorr) return;
 438
 439   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 440     delete fLookUpDxCorr[iPhi];
 441     delete fLookUpDyCorr[iPhi];
 442     delete fLookUpDzCorr[iPhi];
 443
 444     delete fLookUpDxDist[iPhi];
 445     delete fLookUpDyDist[iPhi];
 446     delete fLookUpDzDist[iPhi];
 447   }
 448
 449   delete [] fLookUpDxCorr;
 450   delete [] fLookUpDyCorr;
 451   delete [] fLookUpDzCorr;
 452
 453   delete [] fLookUpDxDist;
 454   delete [] fLookUpDyDist;
 455   delete [] fLookUpDzDist;
 456
 457   fLookUpDxCorr   = 0x0;
 458   fLookUpDyCorr = 0x0;
 459   fLookUpDzCorr    = 0x0;
 460
 461   fLookUpDxDist   = 0x0;
 462   fLookUpDyDist = 0x0;
 463   fLookUpDzDist    = 0x0;
 464 }
 465
 466 //_________________________________________________________________________________________
 467 void AliTPCCorrectionLookupTable::SetupDefaultLimits()
 468 {
 469   //
 470   // Set default limits for tables
 471   //
 472
 473   fNR   = kNR;
 474   fNPhi = kNPhi;
 475   fNZ   = kNZ;
 476   fLimitsR.  ResizeTo(fNR);
 477   fLimitsPhi.ResizeTo(fNPhi);
 478   fLimitsZ.  ResizeTo(fNZ);
 479   fLimitsR.  SetElements(fgkRList);
 480   fLimitsPhi.SetElements(fgkPhiList);
 481   fLimitsZ.  SetElements(fgkZList);
 482 }
 483
 484 //_________________________________________________________________________________________
 485 void AliTPCCorrectionLookupTable::ResetLimits()
 486 {
 487   fNR   = 0;
 488   fNPhi = 0;
 489   fNZ   = 0;
 490
 491   fLimitsR.  ResizeTo(1);
 492   fLimitsPhi.ResizeTo(1);
 493   fLimitsZ.  ResizeTo(1);
 494 }
 495
 496 //_________________________________________________________________________________________
 497 void AliTPCCorrectionLookupTable::MergePhiTables(const char* files)
 498 {
 499   //
 500   // merge all lookup tables stored in 'files' with this one
 501   // assume that each lookup table in each file has only one phi bin
 502   //
 503
 504   ResetTables();
 505   ResetLimits(); // use limits from the first file assuming they are all the same
 506
 507   TString sfiles=gSystem->GetFromPipe(Form("ls %s",files));
 508   TObjArray *arrFiles=sfiles.Tokenize("\n");
 509
 510   for (Int_t i=0; i<arrFiles->GetEntriesFast();++i){
 511     TFile f(arrFiles->At(i)->GetName());
 512     if (!f.IsOpen() || f.IsZombie()) continue;
 513     AliTPCCorrectionLookupTable *lt=dynamic_cast<AliTPCCorrectionLookupTable*>(f.Get(f.GetListOfKeys()->First()->GetName()));
 514     if (!lt) {
 515       AliError(Form("first object in file '%s' is not of type AliTPCCorrectionLookupTable!",f.GetName()));
 516       continue;
 517     }
 518
 519     if (!fNR) {
 520       fNR        = lt->fNR;
 521       fNPhi      = lt->fNPhi;
 522       fNZ        = lt->fNZ;
 523       fLimitsR   = lt->fLimitsR;
 524       fLimitsZ   = lt->fLimitsZ;
 525       fLimitsPhi = lt->fLimitsPhi;
 526       InitTableArrays();
 527     } else {
 528       if (fNR!=lt->fNR || fNPhi!=lt->fNPhi || fNZ!=lt->fNZ ){
 529         AliError(Form("Current limits don't macht the ones in file '%s'",f.GetName()));
 530         continue;
 531       }
 532     }
 533
 534     for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi) {
 535       if (!lt->fLookUpDxCorr[iPhi]) continue;
 536
 537       AliInfo(Form("Adding phi bin '%d' from file '%s'",iPhi,f.GetName()));
 538
 539       InitTablesPhiBin(iPhi);
 540
 541       *fLookUpDxDist[iPhi] = *lt->fLookUpDxDist[iPhi];
 542       *fLookUpDyDist[iPhi] = *lt->fLookUpDyDist[iPhi];
 543       *fLookUpDzDist[iPhi] = *lt->fLookUpDzDist[iPhi];
 544       //
 545       *fLookUpDxCorr[iPhi] = *lt->fLookUpDxCorr[iPhi];
 546       *fLookUpDyCorr[iPhi] = *lt->fLookUpDyCorr[iPhi];
 547       *fLookUpDzCorr[iPhi] = *lt->fLookUpDzCorr[iPhi];
 548       break;
 549     }
 550   }
 551
 552   //check of all phi bins are initialised
 553   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi) {
 554     if (!fLookUpDxCorr[iPhi]) {
 555       AliFatal(Form("Phi bin '%d' not initialised from files!",iPhi));
 556     }
 557   }
 558
 559   delete arrFiles;
 560 }
 561
 562 //_________________________________________________________________________________________
 563 void AliTPCCorrectionLookupTable::BuildExactInverse()
 564 {
 565   //
 566   // this method build the exact inverse of the standard distortion map
 567   // for the the distortion man first needs to be calculated
 568   // then the correction map will be overwritten
 569   //
 570
 571   Float_t x[3]    = {0.,0.,0.};
 572   Float_t x2[3]   = {0.,0.,0.};
 573   Float_t xref[3] = {0.,0.,0.};
 574   Float_t xd[3]   = {0.,0.,0.};
 575   Float_t dx[3]   = {0.,0.,0.};
 576
 577   // reset correction matrices
 578   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 579     TMatrixF &mDxCorr   = *fLookUpDxCorr[iPhi];
 580     TMatrixF &mDyCorr   = *fLookUpDyCorr[iPhi];
 581     TMatrixF &mDzCorr   = *fLookUpDzCorr[iPhi];
 582
 583     for (Int_t ir=0; ir<fNR; ++ir){
 584       for (Int_t iz=0; iz<fNZ; ++iz){
 585         mDxCorr(ir,iz) = -1000.;
 586         mDyCorr(ir,iz) = -1000.;
 587         mDzCorr(ir,iz) = -1000.;
 588       }
 589     }
 590   }
 591
 592   // get interplolated corrections on standard grid
 593   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 594     Double_t phi=fLimitsPhi(iPhi);
 595     TMatrixF &mDxDist   = *fLookUpDxDist[iPhi];
 596     TMatrixF &mDyDist   = *fLookUpDyDist[iPhi];
 597     TMatrixF &mDzDist   = *fLookUpDzDist[iPhi];
 598
 599     for (Int_t ir=0; ir<fNR; ++ir){
 600       Double_t r=fLimitsR(ir);
 601       x[0]=r * TMath::Cos(phi);
 602       x[1]=r * TMath::Sin(phi);
 603
 604       for (Int_t iz=0; iz<fNZ; ++iz){
 605         Double_t z=fLimitsZ(iz);
 606         x[2]=z;
 607
 608         //TODO: change hardcoded value for r>133.?
 609         Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
 610         if (r>133.) roc+=36;
 611         if (z<0)    roc+=18;
 612
 613         dx[0] = mDxDist(ir,iz);
 614         dx[1] = mDyDist(ir,iz);
 615         dx[2] = mDzDist(ir,iz);
 616
 617         xd[0] = x[0]+dx[0];
 618         xd[1] = x[1]+dx[1];
 619         xd[2] = x[2]+dx[2];
 620
 621         const Double_t phid = TVector2::Phi_0_2pi(TMath::ATan2(xd[1],xd[0]));
 622         const Double_t rd   = TMath::Sqrt(xd[0]*xd[0] + xd[1]*xd[1]);
 623         const Double_t zd   = xd[2];
 624
 625         Int_t ilow = 0, jlow = 0, klow = 0 ;
 626
 627         Search( fLimitsR.GetNrows(),   fLimitsR.GetMatrixArray(),   rd,   ilow   ) ;
 628         Search( fLimitsZ.GetNrows(),   fLimitsZ.GetMatrixArray(),   zd,   jlow   ) ;
 629         Search( fLimitsPhi.GetNrows(), fLimitsPhi.GetMatrixArray(), phid, klow   ) ;
 630
 631         if ( ilow < 0 ) ilow = 0 ;   // check if out of range
 632         if ( jlow < 0 ) jlow = 0 ;
 633         if ( klow < 0 ) klow = 0 ;
 634         if ( ilow >= fLimitsR.GetNrows())   ilow = fLimitsR.GetNrows() - 1;
 635         if ( jlow >= fLimitsZ.GetNrows())   jlow = fLimitsZ.GetNrows() - 1;
 636         if ( klow >= fLimitsPhi.GetNrows()) klow = fLimitsPhi.GetNrows() - 1;
 637
 638         const Double_t phiRef = fLimitsPhi[klow];
 639         const Double_t rRef   = fLimitsR[ilow];
 640         const Double_t zRef   = fLimitsZ[jlow];
 641
 642         TMatrixF &mDxCorr   = *fLookUpDxCorr[klow];
 643         if ( mDxCorr(ilow, jlow) > -1000. ) continue;
 644         TMatrixF &mDyCorr   = *fLookUpDyCorr[klow];
 645         TMatrixF &mDzCorr   = *fLookUpDzCorr[klow];
 646
 647         xref[0]= rRef * TMath::Cos(phiRef);
 648         xref[1]= rRef * TMath::Sin(phiRef);
 649         xref[2]= zRef;
 650
 651         FindClosestPosition(ir,iz,iPhi, xref, x2);
 652
 653         GetDistortion(x2,roc,dx);
 654
 655         mDxCorr(ilow, jlow) = -dx[0];
 656         mDyCorr(ilow, jlow) = -dx[1];
 657         mDzCorr(ilow, jlow) = -dx[2];
 658
 659 //         printf("%3d %3d %3d\n",iPhi, ir, iz);
 660 //         printf("%3d %3d %3d\n",klow, ilow, jlow);
 661 //         printf("x2:   %.5f %.5f %.5f\n", x2[0], x2[1], x2[2]);
 662 //         printf("x2d:  %.5f %.5f %.5f\n", x2[0]+dx[0], x2[1]+dx[1], x2[2]+dx[2]);
 663 //         printf("xref: %.5f %.5f %.5f\n", xref[0], xref[1], xref[2]);
 664 //         printf("xrd:  %.5f %.5f %.5f\n", x2[0]+dx[0]-xref[0], x2[1]+dx[1]-xref[1], x2[2]+dx[2]-xref[2]);
 665 //         printf("phid: %.5f %.5f %.5f\n", phid,rd,zd);
 666 //         printf("phir: %.5f %.5f %.5f\n", phiRef,rRef,zRef);
 667 //         printf("\n");
 668       }
 669     }
 670   }
 671
 672   // fill remaining empty bins
 673   // The last ein first phi bin entries must be identical, fill those first
 674   {
 675   TMatrixF &mDxCorr   = *fLookUpDxCorr[0];
 676   TMatrixF &mDyCorr   = *fLookUpDyCorr[0];
 677   TMatrixF &mDzCorr   = *fLookUpDzCorr[0];
 678
 679   TMatrixF &mDxCorr2  = *fLookUpDxCorr[fNPhi-1];
 680   TMatrixF &mDyCorr2  = *fLookUpDyCorr[fNPhi-1];
 681   TMatrixF &mDzCorr2  = *fLookUpDzCorr[fNPhi-1];
 682
 683   for (Int_t ir=0; ir<fNR; ++ir){
 684     for (Int_t iz=0; iz<fNZ; ++iz){
 685       mDxCorr2(ir,iz) = mDxCorr(ir,iz);
 686       mDyCorr2(ir,iz) = mDyCorr(ir,iz);
 687       mDzCorr2(ir,iz) = mDzCorr(ir,iz);
 688     }
 689   }
 690   }
 691
 692   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 693     TMatrixF &mDxCorr   = *fLookUpDxCorr[iPhi];
 694     TMatrixF &mDyCorr   = *fLookUpDyCorr[iPhi];
 695     TMatrixF &mDzCorr   = *fLookUpDzCorr[iPhi];
 696
 697     Double_t phi=fLimitsPhi(iPhi);
 698     for (Int_t ir=0; ir<fNR; ++ir){
 699       Double_t r=fLimitsR(ir);
 700       x[0]=r * TMath::Cos(phi);
 701       x[1]=r * TMath::Sin(phi);
 702
 703       for (Int_t iz=0; iz<fNZ; ++iz){
 704         if (mDxCorr(ir,iz) > -999.) continue;
 705
 706         Double_t z=fLimitsZ(iz);
 707         x[2]=z;
 708
 709         //TODO: change hardcoded value for r>133.?
 710         Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
 711         if (r>133.) roc+=36;
 712         if (z<0)    roc+=18;
 713
 714         // get last point
 715         dx[0] = mDxCorr(ir,iz-1);
 716         dx[1] = mDyCorr(ir,iz-1);
 717         dx[2] = mDzCorr(ir,iz-1);
 718
 719         xd[0] = x[0]+dx[0];
 720         xd[1] = x[1]+dx[1];
 721         xd[2] = x[2]+dx[2];
 722
 723         // get distorted point
 724         const Double_t phid = TVector2::Phi_0_2pi(TMath::ATan2(xd[1],xd[0]));
 725         const Double_t rd   = TMath::Sqrt(xd[0]*xd[0] + xd[1]*xd[1]);
 726         const Double_t zd   = xd[2];
 727
 728         Int_t ilow = 0, jlow = 0, klow = 0 ;
 729
 730         Search( fLimitsR.GetNrows(),   fLimitsR.GetMatrixArray(),   rd,   ilow   ) ;
 731         Search( fLimitsZ.GetNrows(),   fLimitsZ.GetMatrixArray(),   zd,   jlow   ) ;
 732         Search( fLimitsPhi.GetNrows(), fLimitsPhi.GetMatrixArray(), phid, klow   ) ;
 733
 734         if ( ilow < 0 ) ilow = 0 ;   // check if out of range
 735         if ( jlow < 0 ) jlow = 0 ;
 736         if ( klow < 0 ) klow = 0 ;
 737         if ( ilow >= fLimitsR.GetNrows())   ilow = fLimitsR.GetNrows() - 1;
 738         if ( jlow >= fLimitsZ.GetNrows())   jlow = fLimitsZ.GetNrows() - 1;
 739         if ( klow >= fLimitsPhi.GetNrows()) klow = fLimitsPhi.GetNrows() - 1;
 740
 741         FindClosestPosition(ilow,jlow,klow, x, x2);
 742
 743         GetDistortion(x2,roc,dx);
 744
 745         mDxCorr(ir, iz) = -dx[0];
 746         mDyCorr(ir, iz) = -dx[1];
 747         mDzCorr(ir, iz) = -dx[2];
 748       }
 749     }
 750   }
 751
 752 }
 753
 754 //_________________________________________________________________________________________
 755 void AliTPCCorrectionLookupTable::FindClosestPosition(const Int_t binR, const Int_t binZ, const Int_t binPhi,
 756                                                       const Float_t xref[3], Float_t xret[3])
 757 {
 758   //
 759   //
 760   //
 761
 762 //   static TLinearFitter fitx(2,"pol2");
 763 //   static TLinearFitter fity(2,"pol2");
 764 //   static TLinearFitter fitz(2,"pol2");
 765   static TLinearFitter fitx(4,"hyp3");
 766   static TLinearFitter fity(4,"hyp3");
 767   static TLinearFitter fitz(4,"hyp3");
 768   fitx.ClearPoints();
 769   fity.ClearPoints();
 770   fitz.ClearPoints();
 771
 772   const Int_t nPoints=3;
 773   Int_t counter=0;
 774   Int_t rMin=binR;
 775   Int_t zMin=binZ;
 776   Int_t phiMin=binPhi;
 777
 778   counter=nPoints/2;
 779   while (rMin>0 && counter--) --rMin;
 780   counter=nPoints/2;
 781   while (zMin>0 && counter--) --zMin;
 782   counter=nPoints/2;
 783   while (phiMin>0 && counter--) --phiMin;
 784
 785   Int_t rMax   = rMin  +nPoints;
 786   Int_t zMax   = zMin  +nPoints;
 787   Int_t phiMax = phiMin+nPoints;
 788
 789   while (rMax>=fNR) {--rMin; --rMax;}
 790   while (zMax>=fNZ) {--zMin; --zMax;}
 791   while (phiMax>=fNPhi) {--phiMin; --phiMax;}
 792
 793   Float_t  x[3]    = {0.,0.,0.};
 794   Double_t xd[3]   = {0.,0.,0.};
 795   Float_t  dx[3]   = {0.,0.,0.};
 796
 797   for (Int_t iPhi=phiMin; iPhi<phiMax; ++iPhi) {
 798     TMatrixF &mDxDist   = *fLookUpDxDist[iPhi];
 799     TMatrixF &mDyDist   = *fLookUpDyDist[iPhi];
 800     TMatrixF &mDzDist   = *fLookUpDzDist[iPhi];
 801
 802     Double_t phi=fLimitsPhi(iPhi);
 803     for (Int_t ir=rMin; ir<rMax; ++ir){
 804       Double_t r=fLimitsR(ir);
 805       x[0]=r * TMath::Cos(phi);
 806       x[1]=r * TMath::Sin(phi);
 807
 808       for (Int_t iz=zMin; iz<zMax; ++iz){
 809         Double_t z=fLimitsZ(iz);
 810         x[2]=z;
 811
 812         dx[0] = mDxDist(ir,iz);
 813         dx[1] = mDyDist(ir,iz);
 814         dx[2] = mDzDist(ir,iz);
 815
 816         xd[0] = x[0]+dx[0];
 817         xd[1] = x[1]+dx[1];
 818         xd[2] = x[2]+dx[2];
 819
 820         fitx.AddPoint(xd,   x[0]);
 821         fity.AddPoint(xd, x[1]);
 822         fitz.AddPoint(xd, x[2]);
 823       }
 824     }
 825   }
 826
 827   fitx.Eval();
 828   fity.Eval();
 829   fitz.Eval();
 830   xret[0] = fitx.GetParameter(0) + fitx.GetParameter(1)*xref[0]
 831                                  + fitx.GetParameter(2)*xref[1]
 832                                  + fitx.GetParameter(3)*xref[2];
 833   xret[1] = fity.GetParameter(0) + fity.GetParameter(1)*xref[0]
 834                                  + fity.GetParameter(2)*xref[1]
 835                                  + fity.GetParameter(3)*xref[2];
 836   xret[2] = fitz.GetParameter(0) + fitz.GetParameter(1)*xref[0]
 837                                  + fitz.GetParameter(2)*xref[1]
 838                                  + fitz.GetParameter(3)*xref[2];
 839 //   xret[0] = fitx.GetParameter(0) + fitx.GetParameter(1)*xref[0] + fitx.GetParameter(2)*xref[0]*xref[0];
 840 //   xret[1] = fity.GetParameter(0) + fity.GetParameter(1)*xref[1] + fity.GetParameter(2)*xref[1]*xref[1];
 841 //   xret[2] = fitz.GetParameter(0) + fitz.GetParameter(1)*xref[2] + fitz.GetParameter(2)*xref[2]*xref[2];
 842 }
 843