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
  20 #include <AliLog.h>
  21 #include <AliTPCROC.h>
  22
  23 #include "AliTPCCorrection.h"
  24
  25 #include "AliTPCCorrectionLookupTable.h"
  26
  27 ClassImp(AliTPCCorrectionLookupTable)
  28
  29 //_________________________________________________________________________________________
  30 AliTPCCorrectionLookupTable::AliTPCCorrectionLookupTable()
  31 : AliTPCCorrection()
  32 , fNR(90)
  33 , fNPhi(144)
  34 , fNZ(130)
  35 , fLimitsRows()
  36 , fLimitsPhiSlices()
  37 , fLimitsColumns()
  38 , fLookUpDxDist(0x0)
  39 , fLookUpDyDist(0x0)
  40 , fLookUpDzDist(0x0)
  41 , fLookUpDxCorr(0x0)
  42 , fLookUpDyCorr(0x0)
  43 , fLookUpDzCorr(0x0)
  44 {
  45   //
  46   //
  47   //
  48 }
  49 //_________________________________________________________________________________________
  50 AliTPCCorrectionLookupTable::~AliTPCCorrectionLookupTable()
  51 {
  52   //
  53   // dtor
  54   //
  55
  56   ResetTables();
  57 }
  58
  59 //_________________________________________________________________________________________
  60 void AliTPCCorrectionLookupTable::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
  61   //
  62   // Get interpolated Distortion
  63   //
  64
  65   GetInterpolation(x,roc,dx,fLookUpDxDist,fLookUpDyDist,fLookUpDzDist);
  66 }
  67
  68 //_________________________________________________________________________________________
  69 void AliTPCCorrectionLookupTable::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
  70   //
  71   // Get interplolated correction
  72   //
  73   GetInterpolation(x,roc,dx,fLookUpDxCorr,fLookUpDyCorr,fLookUpDzCorr);
  74 }
  75
  76 //_________________________________________________________________________________________
  77 void AliTPCCorrectionLookupTable::GetInterpolation(const Float_t x[],const Short_t roc,Float_t dx[],
  78                                                    TMatrixF **mDx, TMatrixF **mDy, TMatrixF **mDz)
  79 {
  80   //
  81   // Calculates the correction/distotring from a lookup table
  82   // type: 0 = correction
  83   //       1 = distortion
  84   //
  85
  86 //   Float_t typeSign=-1;
  87 //   if (type==1) typeSign=1;
  88
  89   Int_t   order     = 1 ;    // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
  90
  91   Double_t r, phi, z ;
  92   Int_t    sign;
  93
  94   r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
  95   phi    =  TMath::ATan2(x[1],x[0]) ;
  96   if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
  97   z      =  x[2] ;                                         // Create temporary copy of x[2]
  98
  99   if ( (roc%36) < 18 ) {
 100     sign =  1;       // (TPC A side)
 101   } else {
 102     sign = -1;       // (TPC C side)
 103   }
 104
 105   if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
 106   if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
 107
 108
 109   if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
 110     AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
 111
 112   // Get the Er and Ephi field integrals plus the integral over Z
 113     dx[0] = Interpolate3DTable(order, r, z, phi,
 114                                fNR, fNZ, fNPhi,
 115                                fLimitsRows.GetMatrixArray(),
 116                                fLimitsColumns.GetMatrixArray(),
 117                                fLimitsPhiSlices.GetMatrixArray(),
 118                                mDx  );
 119     dx[1] = Interpolate3DTable(order, r, z, phi,
 120                                fNR, fNZ, fNPhi,
 121                                fLimitsRows.GetMatrixArray(),
 122                                fLimitsColumns.GetMatrixArray(),
 123                                fLimitsPhiSlices.GetMatrixArray(),
 124                                mDy);
 125     dx[2] = Interpolate3DTable(order, r, z, phi,
 126                                fNR, fNZ, fNPhi,
 127                                fLimitsRows.GetMatrixArray(),
 128                                fLimitsColumns.GetMatrixArray(),
 129                                fLimitsPhiSlices.GetMatrixArray(),
 130                                mDz   );
 131 }
 132
 133 //_________________________________________________________________________________________
 134 void AliTPCCorrectionLookupTable::CreateLookupTable(AliTPCCorrection &tpcCorr, Float_t stepSize/*=5.*/)
 135 {
 136   //
 137   //
 138   //
 139
 140   // create distortion lookup table
 141
 142   TStopwatch s;
 143
 144   ResetTables();
 145   InitTables();
 146
 147   const Float_t delta=stepSize; // 5cm
 148   Float_t x[3]={0.,0.,0.};
 149   Float_t dx[3]={0.,0.,0.};
 150
 151   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 152     Double_t phi=fLimitsPhiSlices(iPhi);
 153     //
 154     TMatrixF &mDxDist   = *fLookUpDxDist[iPhi];
 155     TMatrixF &mDyDist   = *fLookUpDyDist[iPhi];
 156     TMatrixF &mDzDist   = *fLookUpDzDist[iPhi];
 157     //
 158     TMatrixF &mDxCorr   = *fLookUpDxCorr[iPhi];
 159     TMatrixF &mDyCorr   = *fLookUpDyCorr[iPhi];
 160     TMatrixF &mDzCorr   = *fLookUpDzCorr[iPhi];
 161
 162     for (Int_t ir=0; ir<fNR; ++ir){
 163       Double_t r=fLimitsRows(ir);
 164       x[0]=r * TMath::Cos(phi);
 165       x[1]=r * TMath::Sin(phi);
 166
 167       for (Int_t iz=0; iz<fNZ; ++iz){
 168         Double_t z=fLimitsColumns(iz);
 169         x[2]=z;
 170         //TODO: change hardcoded value for r>133.?
 171         Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
 172         if (r>133.) roc+=36;
 173         if (z<0)    roc+=18;
 174
 175         tpcCorr.GetDistortionIntegralDz(x,roc,dx,delta);
 176         mDxDist(ir,iz)=dx[0];
 177         mDyDist(ir,iz)=dx[1];
 178         mDzDist(ir,iz)=dx[2];
 179
 180         tpcCorr.GetCorrectionIntegralDz(x,roc,dx,delta);
 181         mDxCorr(ir,iz)=dx[0];
 182         mDyCorr(ir,iz)=dx[1];
 183         mDzCorr(ir,iz)=dx[2];
 184
 185       }
 186     }
 187   }
 188
 189   s.Stop();
 190   AliInfo("Required time for lookup table creation: %.2f (%.2f) sec. real (cpu)");
 191 }
 192
 193 //_________________________________________________________________________________________
 194 void AliTPCCorrectionLookupTable::InitTables()
 195 {
 196   //
 197   //
 198   //
 199   fLookUpDxCorr = new TMatrixF*[fNPhi];
 200   fLookUpDyCorr = new TMatrixF*[fNPhi];
 201   fLookUpDzCorr = new TMatrixF*[fNPhi];
 202
 203   fLookUpDxDist = new TMatrixF*[fNPhi];
 204   fLookUpDyDist = new TMatrixF*[fNPhi];
 205   fLookUpDzDist = new TMatrixF*[fNPhi];
 206
 207   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 208     fLookUpDxCorr[iPhi] = new TMatrixF(fNR,fNZ);
 209     fLookUpDyCorr[iPhi] = new TMatrixF(fNR,fNZ);
 210     fLookUpDzCorr[iPhi] = new TMatrixF(fNR,fNZ);
 211
 212     fLookUpDxDist[iPhi] = new TMatrixF(fNR,fNZ);
 213     fLookUpDyDist[iPhi] = new TMatrixF(fNR,fNZ);
 214     fLookUpDzDist[iPhi] = new TMatrixF(fNR,fNZ);
 215   }
 216
 217   SetupDefaultLimits();
 218 }
 219
 220 //_________________________________________________________________________________________
 221 void AliTPCCorrectionLookupTable::ResetTables()
 222 {
 223   //
 224   // Reset the lookup tables
 225   //
 226
 227   if (!fLookUpDxCorr) return;
 228
 229   for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
 230     delete fLookUpDxCorr[iPhi];
 231     delete fLookUpDyCorr[iPhi];
 232     delete fLookUpDzCorr[iPhi];
 233
 234     delete fLookUpDxDist[iPhi];
 235     delete fLookUpDyDist[iPhi];
 236     delete fLookUpDzDist[iPhi];
 237   }
 238
 239   delete [] fLookUpDxCorr;
 240   delete [] fLookUpDyCorr;
 241   delete [] fLookUpDzCorr;
 242
 243   delete [] fLookUpDxDist;
 244   delete [] fLookUpDyDist;
 245   delete [] fLookUpDzDist;
 246
 247   fLookUpDxCorr   = 0x0;
 248   fLookUpDyCorr = 0x0;
 249   fLookUpDzCorr    = 0x0;
 250
 251   fLookUpDxDist   = 0x0;
 252   fLookUpDyDist = 0x0;
 253   fLookUpDzDist    = 0x0;
 254
 255   fLimitsRows.ResizeTo(1);
 256   fLimitsPhiSlices.ResizeTo(1);
 257   fLimitsColumns.ResizeTo(1);
 258 }
 259
 260 //_________________________________________________________________________________________
 261 void AliTPCCorrectionLookupTable::SetupDefaultLimits()
 262 {
 263   //
 264   // Set default limits for tables
 265   //
 266
 267   fLimitsRows.ResizeTo(fNR);
 268   fLimitsPhiSlices.ResizeTo(fNPhi);
 269   fLimitsColumns.ResizeTo(fNZ);
 270
 271   Double_t *limRList   = fLimitsRows.GetMatrixArray();
 272   Double_t *limPhiList = fLimitsPhiSlices.GetMatrixArray();
 273   Double_t *limZList   = fLimitsColumns.GetMatrixArray();
 274
 275   AliTPCROC * roc = AliTPCROC::Instance();
 276   const Double_t rLow =  TMath::Floor(roc->GetPadRowRadii(0,0))-1; // first padRow plus some margin
 277
 278   // fulcrums in R
 279   limRList[0] = rLow;
 280   for (Int_t i = 1; i<fNR; i++) {
 281     limRList[i] = limRList[i-1] + 3.5;     // 3.5 cm spacing
 282     if (limRList[i]<90 ||limRList[i]>245){
 283       limRList[i] = limRList[i-1] + 0.5; // 0.5 cm spacing
 284     } else if (limRList[i]<100 || limRList[i]>235){
 285       limRList[i] = limRList[i-1] + 1.5;  // 1.5 cm spacing
 286     } else if (limRList[i]<120 || limRList[i]>225){
 287       limRList[i] = limRList[i-1] + 2.5;  // 2.5 cm spacing
 288     }
 289   }
 290
 291   // fulcrums in Z
 292   limZList[0] = -249.5;
 293   limZList[fNZ-1] = 249.5;
 294   for (Int_t j = 1; j<fNZ/2; j++) {
 295     limZList[j] = limZList[j-1];
 296     if      (TMath::Abs(limZList[j])< 0.15){
 297       limZList[j] = limZList[j-1] + 0.09; // 0.09 cm spacing
 298     } else if(TMath::Abs(limZList[j])< 0.6){
 299       limZList[j] = limZList[j-1] + 0.4; // 0.4 cm spacing
 300     } else if      (TMath::Abs(limZList[j])< 2.5 || TMath::Abs(limZList[j])>248){
 301       limZList[j] = limZList[j-1] + 0.5; // 0.5 cm spacing
 302     } else if (TMath::Abs(limZList[j])<10 || TMath::Abs(limZList[j])>235){
 303       limZList[j] = limZList[j-1] + 1.5;  // 1.5 cm spacing
 304     } else if (TMath::Abs(limZList[j])<25 || TMath::Abs(limZList[j])>225){
 305       limZList[j] = limZList[j-1] + 2.5;  // 2.5 cm spacing
 306     } else{
 307       limZList[j] = limZList[j-1] + 4;  // 4 cm spacing
 308     }
 309
 310     limZList[fNZ-j-1] = -limZList[j];
 311   }
 312
 313   // fulcrums in phi
 314   for (Int_t k = 0; k<fNPhi; k++)
 315     limPhiList[k] = TMath::TwoPi()*k/(fNPhi-1);
 316
 317 }