TPC/CalibMacros/AliTPCDistortions.cxx

   1 //
   2 //
   3 // This class is based on Jim Thomas's implementation of the ExB effects.
   4 // It was shortened and made to comply with AliROOT coding conventions,
   5 // but not changed in functionallity.
   6 //
   7 //
   8
   9
  10 /***********************************************************************
  11  *
  12  * Author: Jim Thomas   11/4/2009
  13  *
  14  ***********************************************************************
  15  *
  16  * Description:
  17  *    Utilities for handling SpaceCharge and a few other distortions
  18  *    that may occur in the ALICE TPC
  19  *
  20  ***********************************************************************/
  21
  22 #include "AliTPCDistortions.h"
  23 #include <AliMagF.h>
  24 #include <TMath.h>
  25
  26 //
  27 //  Standard maps for E and B Field Distortions
  28 //
  29 //  Note the careful steps in Z around the Central Electrode due to possible discontinuities at CE.
  30 //  Needed for interpolation tools on the grid.  Also note that whenever we interpolate this grid,
  31 //  we explicitly do not allow Z to get closer to CE than 0.2 cm.  This gives three points for
  32 //  quadratic interpolation in all cases (if you need it).
  33
  34 const Double_t AliTPCDistortions::fgkRList[AliTPCDistortions::kNR] =  {   84.0,   84.5,   85.0,   85.5,   86.0,  87.0,    88.0,
  35                                    90.0,   92.0,   94.0,   96.0,   98.0,  100.0,  102.0,  104.0,  106.0,  108.0,
  36                                   110.0,  112.0,  114.0,  116.0,  118.0,  120.0,  122.0,  124.0,  126.0,  128.0,
  37                                   130.0,  132.0,  134.0,  136.0,  138.0,  140.0,  142.0,  144.0,  146.0,  148.0,
  38                                   150.0,  152.0,  154.0,  156.0,  158.0,  160.0,  162.0,  164.0,  166.0,  168.0,
  39                                   170.0,  172.0,  174.0,  176.0,  178.0,  180.0,  182.0,  184.0,  186.0,  188.0,
  40                                   190.0,  192.0,  194.0,  196.0,  198.0,  200.0,  202.0,  204.0,  206.0,  208.0,
  41                                   210.0,  212.0,  214.0,  216.0,  218.0,  220.0,  222.0,  224.0,  226.0,  228.0,
  42                                   230.0,  232.0,  234.0,  236.0,  238.0,  240.0,  242.0,  244.0,  246.0,  248.0,
  43                                   249.0,  249.5,  250.0,  251.5,  252.0  } ;
  44
  45 const Double_t AliTPCDistortions::fgkPhiList[AliTPCDistortions::kNPhi] = {  0.0000,
  46                                   2.0*TMath::Pi()* 1.0/18, 2.0*TMath::Pi()* 2.0/18, 2.0*TMath::Pi()* 3.0/18,
  47                                   2.0*TMath::Pi()* 4.0/18, 2.0*TMath::Pi()* 5.0/18, 2.0*TMath::Pi()* 6.0/18,
  48                                   2.0*TMath::Pi()* 7.0/18, 2.0*TMath::Pi()* 8.0/18, 2.0*TMath::Pi()* 9.0/18,
  49                                   2.0*TMath::Pi()*10.0/18, 2.0*TMath::Pi()*11.0/18, 2.0*TMath::Pi()*12.0/18,
  50                                   2.0*TMath::Pi()*13.0/18, 2.0*TMath::Pi()*14.0/18, 2.0*TMath::Pi()*15.0/18,
  51                                   2.0*TMath::Pi()*16.0/18, 2.0*TMath::Pi()*17.0/18, 2.0*TMath::Pi() } ;
  52                                   // 19 planes of phi - 18+1 so can wrap around
  53
  54 const Double_t AliTPCDistortions::fgkZList[AliTPCDistortions::kNZ]     =
  55                                { -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
  56                                  -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
  57                                  -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
  58                                  -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
  59                                  -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
  60                                  -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
  61                                  -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
  62                                  -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
  63                                  -100.0,  -98.0,  -96.0,  -94.0,  -92.0,  -90.0,  -88.0,  -86.0,  -84.0,  -82.0,
  64                                   -80.0,  -78.0,  -76.0,  -74.0,  -72.0,  -70.0,  -68.0,  -66.0,  -64.0,  -62.0,
  65                                   -60.0,  -58.0,  -56.0,  -54.0,  -52.0,  -50.0,  -48.0,  -46.0,  -44.0,  -42.0,
  66                                   -40.0,  -38.0,  -36.0,  -34.0,  -32.0,  -30.0,  -28.0,  -26.0,  -24.0,  -22.0,
  67                                   -20.0,  -18.0,  -16.0,  -14.0,  -12.0,  -10.0,   -8.0,   -6.0,   -4.0,   -2.0,
  68                                    -1.0,   -0.5,   -0.2,   -0.1,  -0.05,   0.05,    0.1,    0.2,    0.5,    1.0,
  69                                     2.0,    4.0,    6.0,    8.0,   10.0,   12.0,   14.0,   16.0,   18.0,   20.0,
  70                                    22.0,   24.0,   26.0,   28.0,   30.0,   32.0,   34.0,   36.0,   38.0,   40.0,
  71                                    42.0,   44.0,   46.0,   48.0,   50.0,   52.0,   54.0,   56.0,   58.0,   60.0,
  72                                    62.0,   64.0,   66.0,   68.0,   70.0,   72.0,   74.0,   76.0,   78.0,   80.0,
  73                                    82.0,   84.0,   86.0,   88.0,   90.0,   92.0,   94.0,   96.0,   98.0,  100.0,
  74                                   102.0,  104.0,  106.0,  108.0,  110.0,  112.0,  114.0,  116.0,  118.0,  120.0,
  75                                   122.0,  124.0,  126.0,  128.0,  130.0,  132.0,  134.0,  136.0,  138.0,  140.0,
  76                                   142.0,  144.0,  146.0,  148.0,  150.0,  152.0,  154.0,  156.0,  158.0,  160.0,
  77                                   162.0,  164.0,  166.0,  168.0,  170.0,  172.0,  174.0,  176.0,  178.0,  180.0,
  78                                   182.0,  184.0,  186.0,  188.0,  190.0,  192.0,  194.0,  196.0,  198.0,  200.0,
  79                                   202.0,  204.0,  206.0,  208.0,  210.0,  212.0,  214.0,  216.0,  218.0,  220.0,
  80                                   222.0,  224.0,  226.0,  228.0,  230.0,  232.0,  234.0,  236.0,  238.0,  240.0,
  81                                   242.0,  243.0,  244.0,  245.0,  246.0,  247.0,  248.0,  248.5,  249.0,  249.5   } ;
  82
  83 const Double_t AliTPCDistortions::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min,  83.70 max) (cm)
  84 const Double_t AliTPCDistortions::fgkOFCRadius=254.5;  // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
  85 const Double_t AliTPCDistortions::fgkTPC_Z0=249.7;     // Z location of TPC Gated Grid (cm)
  86 const Double_t AliTPCDistortions::fgkZOffSet= 0.2;     // Offset from CE: calculate all distortions closer to CE as if at this point
  87 const Double_t AliTPCDistortions::fgkCathodeV   =-100000.0; // Cathode Voltage (volts)
  88 const Double_t AliTPCDistortions::fgkGG         =-70.0;     // Gating Grid voltage (volts)
  89 const Double_t AliTPCDistortions::fgkAliceDriftV= 2.73;     // Drift Velocity (cm/microSec) Magnitude
  90
  91
  92
  93
  94 AliTPCDistortions* AliTPCDistortions::fgInstance = 0;
  95
  96
  97
  98
  99 //_ singleton implementation __________________________________________________
 100 AliTPCDistortions* AliTPCDistortions::Instance()
 101 {
 102   //
 103   // Singleton implementation
 104   // Returns an instance of this class, it is created if neccessary
 105   //
 106   if (fgInstance == 0){
 107     fgInstance = new AliTPCDistortions();
 108     //fgInstance->Init();
 109   }
 110   return fgInstance;
 111 }
 112
 113
 114 AliTPCDistortions::AliTPCDistortions()
 115   :fOmegaTau(0.),fT1(1.),fT2(1.),fC0(1.),fC1(0.),fC2(0.),
 116    fBField(0),fXTwist(0.),fYTwist(0.),fIFCShift(0.),fDeltaVGGA(0.),fDeltaVGGC(0.),
 117    fJLow(0),fKLow(0)
 118 {
 119   //
 120   // default constructor (all corrections are set to have zero effect)
 121   //
 122
 123   InitIFCShiftDistortion();
 124   InitGGVoltErrorDistortion();
 125 }
 126
 127 void AliTPCDistortions::RecalculateCs() {
 128   //
 129   // helper function to recalculate the coefficients C0-C2
 130   // after a change of omega-tau or the tensor terms
 131   //
 132   fC0=1.                          / ( 1. + fT2*fT2*fOmegaTau*fOmegaTau ) ;
 133   fC1=fT1*fOmegaTau               / ( 1. + fT1*fT1*fOmegaTau*fOmegaTau ) ;
 134   fC2=fT2*fT2*fOmegaTau*fOmegaTau / ( 1. + fT2*fT2*fOmegaTau*fOmegaTau ) ;
 135 }
 136
 137 void AliTPCDistortions::UndoTwistDistortion(const Double_t x[],Double_t xprime[],Int_t roc)
 138 {
 139   //
 140   // Twist distortion
 141   //
 142   // Remove the effects of a simple "twist" of the TPC in the magnet.  If there is
 143   // an angle between the E and B fields, there will be a distortion in the recorded
 144   // tracks.  This routine takes out that distortion.
 145   //
 146   // The parameters are:
 147   // twist[0]: twist angle in x-z
 148   // twist[1]: twist angle in y-z
 149   //
 150
 151   Double_t        zdrift ;
 152   Int_t           sign ;
 153
 154   Double_t z = x[2] ;                                 // Creat temporary copy of x[2]
 155
 156   if ( roc <  36 ) sign =  1 ;                       // (TPC End A)
 157   else             sign = -1 ;                       // (TPC End C)
 158   if ( roc <= 35 && z <  fgkZOffSet ) z =  fgkZOffSet ;    // Protect against discontinuity at CE
 159   if ( roc >= 36 && z > -fgkZOffSet ) z = -fgkZOffSet ;    // Protect against discontinuity at CE
 160
 161   zdrift    = sign * ( fgkTPC_Z0 - TMath::Abs(z) ) ;
 162
 163   xprime[0] = x[0] - (     fC1 * fYTwist - fC2 * fXTwist ) * zdrift/1000 ;
 164   xprime[1] = x[1] - ( -1* fC1 * fXTwist - fC2 * fYTwist ) * zdrift/1000 ;
 165   xprime[2] = x[2] ;                                       // Subtract to undo the distortion
 166 }
 167
 168
 169 void AliTPCDistortions::InitIFCShiftDistortion() {
 170   //      cout << "AliTPCDistortions::IFCShift   Please wait for the tables to fill ...  ~5 seconds" << endl ;
 171   Int_t nterms = 100 ;
 172   Double_t r,z;
 173   for ( Int_t i = 0 ; i < kNZ ; ++i ) {
 174     z = TMath::Abs( fgkZList[i] ) ;
 175     for ( Int_t j = 0 ; j < kNR ; ++j ) {
 176       r = fgkRList[j] ;
 177       fShiftER[i][j] = 0.0 ;
 178       Double_t intz = 0.0 ;
 179       for ( Int_t n = 1 ; n < nterms ; ++n ) {
 180         Double_t k  = (2*n-1) * TMath::Pi() / fgkTPC_Z0 ;
 181         Double_t Cn = -4.0 / ( k * fgkTPC_Z0 ) ;
 182         Double_t numerator =
 183           TMath::BesselK0( k*fgkOFCRadius ) * TMath::BesselI1( k*r ) +
 184           TMath::BesselK1( k*r )         * TMath::BesselI0( k*fgkOFCRadius ) ;
 185         Double_t denominator =
 186           TMath::BesselK0( k*fgkOFCRadius ) * TMath::BesselI0( k*fgkIFCRadius ) -
 187           TMath::BesselK0( k*fgkIFCRadius ) * TMath::BesselI0( k*fgkOFCRadius ) ;
 188         Double_t zterm = 1 + TMath::Cos( k*z ) ;
 189         Double_t qwe = numerator / denominator ;
 190         intz += Cn * zterm * qwe ;
 191         if ( n>10 && fabs(intz)*1.e-10 > fabs(qwe) ) break;
 192       }
 193       if  ( fgkZList[i] < 0 )  intz = -1 * intz ;  // Force AntiSymmetry of solutions in Z
 194       fShiftER[i][j] = intz ;
 195     }
 196   }
 197 }
 198
 199 void AliTPCDistortions::UndoIFCShiftDistortion(const Double_t x[],Double_t xprime[],Int_t roc)
 200 {
 201   //
 202   // IFC Shift Distortion
 203   //
 204   // The Inner field cage of the TPC may not be perfectly aligned with the outer field cage
 205   // of the TPC.  They can be shifted along the Z axis by up to a 1 mm.  This causes a tilting
 206   // of the equi-potential lines inside the TPC and therefore a DCA error at the vertex.
 207   // The distortion is anti-symmetric in Z.
 208   //
 209
 210   Int_t   order     = 1     ;               // Linear interpolation = 1, Quadratic = 2
 211
 212   Double_t  intEr, intEphi ;
 213   Double_t r, phi, z ;
 214   Int_t    sign ;
 215
 216   r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
 217   phi    =  TMath::ATan2(x[1],x[0]) ;
 218   if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
 219   z      =  x[2] ;                                         // Create temporary copy of x[2]
 220
 221   if ( roc <  36 ) sign =  1 ;                       // (TPC End A)
 222   else                   sign = -1 ;                       // (TPC End C)
 223   if ( roc <= 35 && z <  fgkZOffSet ) z =  fgkZOffSet ;    // Protect against discontinuity at CE
 224   if ( roc >= 36 && z > -fgkZOffSet ) z = -fgkZOffSet ;    // Protect against discontinuity at CE
 225
 226   Interpolate2DEdistortion( order, r, z, fShiftER, intEr ) ;
 227   intEphi = 0.0 ;  // Efield is symmetric in phi
 228
 229   // Subtract to Undo the distortions
 230   if ( r > 0.0 )
 231     {
 232       phi =  phi - fIFCShift*( fC0*intEphi - fC1*intEr ) / r ;
 233       r   =  r   - fIFCShift*( fC0*intEr   + fC1*intEphi ) ;
 234     }
 235
 236   xprime[0] = r * TMath::Cos(phi) ;
 237   xprime[1] = r * TMath::Sin(phi) ;
 238   xprime[2] = x[2] ;
 239
 240 }
 241
 242 void AliTPCDistortions::InitGGVoltErrorDistortion() {
 243   //      cout << "AliTPCDistortions::UndoGG VE  Please wait for the tables to fill ...  ~5 seconds" << endl ;
 244   Double_t r,z;
 245   Int_t nterms = 100 ;
 246   for ( Int_t i = 0 ; i < kNZ ; ++i ) {
 247     z = fgkZList[i] ;
 248     for ( Int_t j = 0 ; j < kNR ; ++j ) {
 249       r = fgkRList[j] ;
 250       fGGVoltErrorER[i][j] = 0.0 ;
 251       Double_t intz = 0.0 ;
 252       for ( Int_t n = 1 ; n < nterms ; ++n ) {
 253         Double_t k    =  n * TMath::Pi() / fgkTPC_Z0 ;
 254         Double_t ein  =  0 ;                    // Error potential on the IFC
 255         Double_t eout =  0 ;                    // Error potential on the OFC
 256         if ( z < 0 ) {
 257           ein   =  2.0 * ((z-fgkTPC_Z0)/fgkTPC_Z0) / ( k * (fgkCathodeV - fgkGG) ) ;
 258           eout  =  2.0 * ((z-fgkTPC_Z0)/fgkTPC_Z0) / ( k * (fgkCathodeV - fgkGG) ) ;
 259         }
 260         if ( z == 0 ) continue ;
 261         if ( z > 0 ) {
 262           ein   =  2.0 * ((fgkTPC_Z0-z)/fgkTPC_Z0) / ( k * (fgkCathodeV - fgkGG) ) ;
 263           eout  =  2.0 * ((fgkTPC_Z0-z)/fgkTPC_Z0) / ( k * (fgkCathodeV - fgkGG) ) ;
 264         }
 265         Double_t an   =  ein  * TMath::BesselK0( k*fgkOFCRadius ) - eout * TMath::BesselK0( k*fgkIFCRadius ) ;
 266         Double_t bn   =  eout * TMath::BesselI0( k*fgkIFCRadius ) - ein  * TMath::BesselI0( k*fgkOFCRadius ) ;
 267         Double_t numerator =
 268           an * TMath::BesselI1( k*r ) - bn * TMath::BesselK1( k*r ) ;
 269         Double_t denominator =
 270           TMath::BesselK0( k*fgkOFCRadius ) * TMath::BesselI0( k*fgkIFCRadius ) -
 271           TMath::BesselK0( k*fgkIFCRadius ) * TMath::BesselI0( k*fgkOFCRadius ) ;
 272         Double_t zterm = TMath::Cos( k*(fgkTPC_Z0-TMath::Abs(z)) ) - 1 ;
 273         intz += zterm * numerator / denominator ;
 274         // Assume series converges, break if small terms
 275         if ( n>10 && fabs(intz)*1.e-10 > fabs(numerator/denominator) ) break;
 276       }
 277       fGGVoltErrorER[i][j] = intz ;
 278     }
 279   }
 280 }
 281
 282
 283
 284 void AliTPCDistortions::UndoGGVoltErrorDistortion(const Double_t x[], Double_t xprime[],Int_t roc)
 285 {
 286   //
 287   // Gated Grid Voltage Error
 288   //
 289   // Calculate the effect of having an incorrect voltage on the A or C end plate Gated Grids.
 290   //
 291   //  Electrostatic Equations from StarNote SN0253 by Howard Wieman.
 292   //  Note that we use a funny coordinate system where Z==0 at the GG.
 293   //
 294
 295   Int_t   order     = 1     ;               // Linear interpolation = 1, Quadratic = 2
 296
 297   Double_t  intEr, intEphi ;
 298   Double_t r, phi, z ;
 299   Int_t    sign ;
 300
 301   Double_t deltaVGG;
 302
 303   r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
 304   phi    =  TMath::ATan2(x[1],x[0]) ;
 305   if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
 306   z      =  x[2] ;
 307
 308   if ( roc <  36 ) sign =  1 ;                       // (TPC End A)
 309   else             sign = -1 ;                       // (TPC End C)
 310   if ( roc <= 35 && z <  fgkZOffSet ) z =  fgkZOffSet ;    // Protect against discontinuity at CE
 311   if ( roc >= 36 && z > -fgkZOffSet ) z = -fgkZOffSet ;    // Protect against discontinuity at CE
 312
 313   if ( roc <  36 ) deltaVGG = fDeltaVGGA ;            // (TPC End A)
 314   else             deltaVGG = fDeltaVGGC ;            // (TPC End C)
 315
 316   Interpolate2DEdistortion( order, r, z, fGGVoltErrorER, intEr ) ;
 317   intEphi = 0.0 ;  // Efield is symmetric in phi
 318
 319   // Subtract to Undo the distortions
 320   if ( r > 0.0 )
 321     {
 322       phi =  phi - deltaVGG*( fC0*intEphi - fC1*intEr ) / r ;
 323       r   =  r   - deltaVGG*( fC0*intEr   + fC1*intEphi ) ;
 324     }
 325
 326   xprime[0] = r * TMath::Cos(phi) ;
 327   xprime[1] = r * TMath::Sin(phi) ;
 328   xprime[2] = x[2] ;
 329
 330 }
 331
 332 void AliTPCDistortions::UndoExBShapeDistortion( const Double_t x[], Double_t xprime[], Int_t roc)
 333 {
 334   //
 335   // Example of using Ruben's integrals to do ExB Shape distortions
 336   //
 337   //
 338
 339   xprime[0] = x[0] ;  xprime[1] = x[1] ;  xprime[2] = x[2]  ;
 340
 341   if (!fBField) return;
 342
 343   Double_t intBStart[3], intBEnd[3] , xStart[3], xEnd[3] ;
 344   Double_t intBxOverBz, intByOverBz;
 345   //Double_t denominator ;
 346   Int_t    sign;
 347   if ( roc <  36 ) sign =  1 ;                    // (TPC End A)
 348   else             sign = -1 ;                    // (TPC End C)
 349
 350   xStart[0] = x[0] ;  xStart[1] = x[1] ;  xStart[2] = x[2]          ;
 351   xEnd  [0] = x[0] ;  xEnd  [1] = x[1] ;  xEnd  [2] = sign * fgkTPC_Z0 ;
 352
 353   fBField -> GetTPCInt(xStart, intBStart) ;
 354   fBField -> GetTPCInt(xEnd  , intBEnd  ) ;
 355
 356   if ( TMath::Abs(intBStart[2] - intBEnd[2]) < 0.1 ) return ;     // Protect against divide by zero, below
 357
 358   intBxOverBz = (intBStart[0] - intBEnd[0]) ;
 359   intByOverBz = (intBStart[1] - intBEnd[1]) ;
 360
 361   xprime[0] +=  ( fC2*intBxOverBz - fC1*intByOverBz )  ;
 362   xprime[1] +=  ( fC2*intByOverBz + fC1*intBxOverBz )  ;
 363   xprime[2] +=  0.0 ;
 364
 365 }
 366
 367
 368 void AliTPCDistortions::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
 369                                                   const Double_t er[kNZ][kNR], Double_t &er_value )
 370 {
 371   //
 372   // Interpolate table - 2D interpolation
 373   //
 374   Double_t save_er[10] ;
 375
 376   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
 377   Search( kNR,   fgkRList,  r,   fKLow   ) ;
 378   if ( fJLow < 0 ) fJLow = 0 ;   // check if out of range
 379   if ( fKLow < 0 ) fKLow = 0 ;
 380   if ( fJLow + order  >=    kNZ - 1 ) fJLow =   kNZ - 1 - order ;
 381   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 382
 383   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ )
 384     {
 385       save_er[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
 386     }
 387   er_value = Interpolate( &fgkZList[fJLow], save_er, order, z )   ;
 388
 389 }
 390
 391
 392 Double_t AliTPCDistortions::Interpolate( const Double_t xArray[], const Double_t yArray[],
 393                                        const Int_t order, const Double_t x )
 394 {
 395   //
 396   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
 397   //
 398
 399   Double_t y ;
 400
 401   if ( order == 2 )                // Quadratic Interpolation = 2
 402
 403     {
 404       y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
 405       y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
 406       y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
 407
 408     }
 409
 410   else                             // Linear Interpolation = 1
 411
 412     {
 413       y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
 414     }
 415
 416   return (y) ;
 417
 418 }
 419
 420
 421 void AliTPCDistortions::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low )
 422 {
 423   //
 424   // Search an ordered table by starting at the most recently used point
 425   //
 426
 427   Long_t middle, high ;
 428   Int_t  ascend = 0, increment = 1 ;
 429
 430   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
 431
 432   if ( low < 0 || low > n-1 ) { low = -1 ; high = n ; }
 433
 434   else                                            // Ordered Search phase
 435     {
 436       if ( (Int_t)( x >= xArray[low] ) == ascend )
 437         {
 438           if ( low == n-1 ) return ;
 439           high = low + 1 ;
 440           while ( (Int_t)( x >= xArray[high] ) == ascend )
 441             {
 442               low = high ;
 443               increment *= 2 ;
 444               high = low + increment ;
 445               if ( high > n-1 )  {  high = n ; break ;  }
 446             }
 447         }
 448       else
 449         {
 450           if ( low == 0 )  {  low = -1 ;  return ;  }
 451           high = low - 1 ;
 452           while ( (Int_t)( x < xArray[low] ) == ascend )
 453             {
 454               high = low ;
 455               increment *= 2 ;
 456               if ( increment >= high )  {  low = -1 ;  break ;  }
 457               else  low = high - increment ;
 458             }
 459         }
 460     }
 461
 462   while ( (high-low) != 1 )                      // Binary Search Phase
 463     {
 464       middle = ( high + low ) / 2 ;
 465       if ( (Int_t)( x >= xArray[middle] ) == ascend )
 466         low = middle ;
 467       else
 468         high = middle ;
 469     }
 470
 471   if ( x == xArray[n-1] ) low = n-2 ;
 472   if ( x == xArray[0]   ) low = 0 ;
 473
 474 }
 475
 476 ClassImp(AliTPCDistortions);