TPC/AliTPCCorrection.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 ////////////////////////////////////////////////////////////////////////////////
  17 //                                                                            //
  18 // AliTPCCorrection class                                                     //
  19 //                                                                            //
  20 // This class provides a general framework to deal with space point           //
  21 // distortions. An correction class which inherits from here is for example   //
  22 // AliTPCExBBShape or AliTPCExBTwist                                          //
  23 //                                                                            //
  24 // General functions are (for example):                                       //
  25 //   CorrectPoint(x,roc) where x is the vector of inital positions in         //
  26 //   cartesian coordinates and roc represents the Read Out chamber number     //
  27 //   according to the offline naming convention. The vector x is overwritten  //
  28 //   with the corrected coordinates.                                          //
  29 //                                                                            //
  30 // An alternative usage would be CorrectPoint(x,roc,dx), which leaves the     //
  31 //   vector x untouched, put returns the distortions via the vector dx        //
  32 //                                                                            //
  33 // The class allows "effective Omega Tau" corrections to be shifted to the    //
  34 // single distortion classes.                                                 //
  35 //                                                                            //
  36 // Note: This class is normally used via the class AliTPCComposedCorrection   //
  37 //                                                                            //
  38 // date: 27/04/2010                                                           //
  39 // Authors: Magnus Mager, Stefan Rossegger, Jim Thomas                        //
  40 ////////////////////////////////////////////////////////////////////////////////
  41 #include "Riostream.h"
  42
  43 #include <TH2F.h>
  44 #include <TMath.h>
  45 #include <TROOT.h>
  46 #include <TTreeStream.h>
  47 #include <TTree.h>
  48 #include <TFile.h>
  49 #include <TTimeStamp.h>
  50 #include <AliCDBStorage.h>
  51 #include <AliCDBId.h>
  52 #include <AliCDBMetaData.h>
  53 #include  "TVectorD.h"
  54
  55
  56 #include "TRandom.h"
  57 #include "AliExternalTrackParam.h"
  58 #include "AliTrackPointArray.h"
  59 #include "TDatabasePDG.h"
  60 #include "AliTrackerBase.h"
  61 #include "AliTPCROC.h"
  62 #include "THnSparse.h"
  63
  64 #include "TRandom.h"
  65 #include "AliTPCTransform.h"
  66 #include "AliTPCcalibDB.h"
  67 #include "AliTPCExB.h"
  68 #include "AliTPCCorrection.h"
  69 #include "AliTPCRecoParam.h"
  70
  71 #include  "AliExternalTrackParam.h"
  72 #include  "AliTrackPointArray.h"
  73 #include  "TDatabasePDG.h"
  74 #include  "AliTrackerBase.h"
  75 #include  "AliTPCROC.h"
  76 #include  "THnSparse.h"
  77
  78 #include  "AliTPCLaserTrack.h"
  79 #include  "AliESDVertex.h"
  80 #include  "AliVertexerTracks.h"
  81 #include  "TDatabasePDG.h"
  82 #include  "TF1.h"
  83
  84 #include "AliTPCCorrection.h"
  85 #include "AliLog.h"
  86
  87 ClassImp(AliTPCCorrection)
  88
  89
  90 TObjArray *AliTPCCorrection::fgVisualCorrection=0;
  91 // instance of correction for visualization
  92
  93
  94 // FIXME: the following values should come from the database
  95 const Double_t AliTPCCorrection::fgkTPCZ0    =249.7;     // nominal gating grid position
  96 const Double_t AliTPCCorrection::fgkIFCRadius= 83.06;    // Mean Radius of the Inner Field Cage ( 82.43 min,  83.70 max) (cm)
  97 const Double_t AliTPCCorrection::fgkOFCRadius=254.5;     // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
  98 const Double_t AliTPCCorrection::fgkZOffSet  = 0.2;      // Offset from CE: calculate all distortions closer to CE as if at this point
  99 const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
 100 const Double_t AliTPCCorrection::fgkGG       =-70.0;     // Gating Grid voltage (volts)
 101
 102
 103
 104 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
 105 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] =  {
 106 84.0,   84.5,   85.0,   85.5,   86.0,  87.0,    88.0,
 107 90.0,   92.0,   94.0,   96.0,   98.0,  100.0,  102.0,  104.0,  106.0,  108.0,
 108 110.0,  112.0,  114.0,  116.0,  118.0,  120.0,  122.0,  124.0,  126.0,  128.0,
 109 130.0,  132.0,  134.0,  136.0,  138.0,  140.0,  142.0,  144.0,  146.0,  148.0,
 110 150.0,  152.0,  154.0,  156.0,  158.0,  160.0,  162.0,  164.0,  166.0,  168.0,
 111 170.0,  172.0,  174.0,  176.0,  178.0,  180.0,  182.0,  184.0,  186.0,  188.0,
 112 190.0,  192.0,  194.0,  196.0,  198.0,  200.0,  202.0,  204.0,  206.0,  208.0,
 113 210.0,  212.0,  214.0,  216.0,  218.0,  220.0,  222.0,  224.0,  226.0,  228.0,
 114 230.0,  232.0,  234.0,  236.0,  238.0,  240.0,  242.0,  244.0,  246.0,  248.0,
 115 249.0,  249.5,  250.0,  251.5,  252.0  } ;
 116
 117 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ]     =   {
 118 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
 119 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
 120 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
 121 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
 122 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
 123 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
 124 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
 125 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
 126 -100.0,  -98.0,  -96.0,  -94.0,  -92.0,  -90.0,  -88.0,  -86.0,  -84.0,  -82.0,
 127 -80.0,  -78.0,  -76.0,  -74.0,  -72.0,  -70.0,  -68.0,  -66.0,  -64.0,  -62.0,
 128 -60.0,  -58.0,  -56.0,  -54.0,  -52.0,  -50.0,  -48.0,  -46.0,  -44.0,  -42.0,
 129 -40.0,  -38.0,  -36.0,  -34.0,  -32.0,  -30.0,  -28.0,  -26.0,  -24.0,  -22.0,
 130 -20.0,  -18.0,  -16.0,  -14.0,  -12.0,  -10.0,   -8.0,   -6.0,   -4.0,   -2.0,
 131 -1.0,   -0.5,   -0.2,   -0.1,  -0.05,   0.05,    0.1,    0.2,    0.5,    1.0,
 132  2.0,    4.0,    6.0,    8.0,   10.0,   12.0,   14.0,   16.0,   18.0,   20.0,
 133  22.0,   24.0,   26.0,   28.0,   30.0,   32.0,   34.0,   36.0,   38.0,   40.0,
 134  42.0,   44.0,   46.0,   48.0,   50.0,   52.0,   54.0,   56.0,   58.0,   60.0,
 135  62.0,   64.0,   66.0,   68.0,   70.0,   72.0,   74.0,   76.0,   78.0,   80.0,
 136  82.0,   84.0,   86.0,   88.0,   90.0,   92.0,   94.0,   96.0,   98.0,  100.0,
 137 102.0,  104.0,  106.0,  108.0,  110.0,  112.0,  114.0,  116.0,  118.0,  120.0,
 138 122.0,  124.0,  126.0,  128.0,  130.0,  132.0,  134.0,  136.0,  138.0,  140.0,
 139 142.0,  144.0,  146.0,  148.0,  150.0,  152.0,  154.0,  156.0,  158.0,  160.0,
 140 162.0,  164.0,  166.0,  168.0,  170.0,  172.0,  174.0,  176.0,  178.0,  180.0,
 141 182.0,  184.0,  186.0,  188.0,  190.0,  192.0,  194.0,  196.0,  198.0,  200.0,
 142 202.0,  204.0,  206.0,  208.0,  210.0,  212.0,  214.0,  216.0,  218.0,  220.0,
 143 222.0,  224.0,  226.0,  228.0,  230.0,  232.0,  234.0,  236.0,  238.0,  240.0,
 144 242.0,  243.0,  244.0,  245.0,  246.0,  247.0,  248.0,  248.5,  249.0,  249.5   } ;
 145
 146
 147
 148 AliTPCCorrection::AliTPCCorrection()
 149   : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
 150 {
 151   //
 152   // default constructor
 153   //
 154   if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
 155 }
 156
 157 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
 158 : TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
 159 {
 160   //
 161   // default constructor, that set the name and title
 162   //
 163   if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
 164 }
 165
 166 AliTPCCorrection::~AliTPCCorrection() {
 167   //
 168   // virtual destructor
 169   //
 170 }
 171
 172 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
 173   //
 174   // Corrects the initial coordinates x (cartesian coordinates)
 175   // according to the given effect (inherited classes)
 176   // roc represents the TPC read out chamber (offline numbering convention)
 177   //
 178   Float_t dx[3];
 179   GetCorrection(x,roc,dx);
 180   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 181 }
 182
 183 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 184   //
 185   // Corrects the initial coordinates x (cartesian coordinates) and stores the new
 186   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 187   // roc represents the TPC read out chamber (offline numbering convention)
 188   //
 189   Float_t dx[3];
 190   GetCorrection(x,roc,dx);
 191   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 192 }
 193
 194 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
 195   //
 196   // Distorts the initial coordinates x (cartesian coordinates)
 197   // according to the given effect (inherited classes)
 198   // roc represents the TPC read out chamber (offline numbering convention)
 199   //
 200   Float_t dx[3];
 201   GetDistortion(x,roc,dx);
 202   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 203 }
 204
 205 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 206   //
 207   // Distorts the initial coordinates x (cartesian coordinates) and stores the new
 208   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 209   // roc represents the TPC read out chamber (offline numbering convention)
 210   //
 211   Float_t dx[3];
 212   GetDistortion(x,roc,dx);
 213   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 214 }
 215
 216 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
 217   //
 218   // This function delivers the correction values dx in respect to the inital coordinates x
 219   // roc represents the TPC read out chamber (offline numbering convention)
 220   // Note: The dx is overwritten by the inherited effectice class ...
 221   //
 222   for (Int_t j=0;j<3;++j) { dx[j]=0.; }
 223 }
 224
 225 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
 226   //
 227   // This function delivers the distortion values dx in respect to the inital coordinates x
 228   // roc represents the TPC read out chamber (offline numbering convention)
 229   //
 230   GetCorrection(x,roc,dx);
 231   for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
 232 }
 233
 234 void AliTPCCorrection::Init() {
 235   //
 236   // Initialization funtion (not used at the moment)
 237   //
 238 }
 239
 240 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
 241   //
 242   // Update function
 243   //
 244 }
 245
 246 void AliTPCCorrection::Print(Option_t* /*option*/) const {
 247   //
 248   // Print function to check which correction classes are used
 249   // option=="d" prints details regarding the setted magnitude
 250   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 251   //
 252   printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
 253 }
 254
 255 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
 256   //
 257   // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
 258   // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
 259   // calibration run
 260   //
 261   fT1=t1;
 262   fT2=t2;
 263   //SetOmegaTauT1T2(omegaTau, t1, t2);
 264 }
 265
 266 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
 267   //
 268   // Simple plot functionality.
 269   // Returns a 2d hisogram which represents the corrections in radial direction (dr)
 270   // in respect to position z within the XY plane.
 271   // The histogramm has nx times ny entries.
 272   //
 273
 274   TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
 275                      nx,-250.,250.,ny,-250.,250.);
 276   Float_t x[3],dx[3];
 277   x[2]=z;
 278   Int_t roc=z>0.?0:18; // FIXME
 279   for (Int_t iy=1;iy<=ny;++iy) {
 280     x[1]=h->GetYaxis()->GetBinCenter(iy);
 281     for (Int_t ix=1;ix<=nx;++ix) {
 282       x[0]=h->GetXaxis()->GetBinCenter(ix);
 283       GetCorrection(x,roc,dx);
 284       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 285       if (90.<=r0 && r0<=250.) {
 286         Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 287         h->SetBinContent(ix,iy,r1-r0);
 288       }
 289       else
 290         h->SetBinContent(ix,iy,0.);
 291     }
 292   }
 293   return h;
 294 }
 295
 296 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
 297   //
 298   // Simple plot functionality.
 299   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 300   // in respect to position z within the XY plane.
 301   // The histogramm has nx times ny entries.
 302   //
 303
 304   TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
 305                      nx,-250.,250.,ny,-250.,250.);
 306   Float_t x[3],dx[3];
 307   x[2]=z;
 308   Int_t roc=z>0.?0:18; // FIXME
 309   for (Int_t iy=1;iy<=ny;++iy) {
 310     x[1]=h->GetYaxis()->GetBinCenter(iy);
 311     for (Int_t ix=1;ix<=nx;++ix) {
 312       x[0]=h->GetXaxis()->GetBinCenter(ix);
 313       GetCorrection(x,roc,dx);
 314       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 315       if (90.<=r0 && r0<=250.) {
 316         Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 317         Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 318
 319         Float_t dphi=phi1-phi0;
 320         if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 321         if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 322
 323         h->SetBinContent(ix,iy,r0*dphi);
 324       }
 325       else
 326         h->SetBinContent(ix,iy,0.);
 327     }
 328   }
 329   return h;
 330 }
 331
 332 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
 333   //
 334   // Simple plot functionality.
 335   // Returns a 2d hisogram which represents the corrections in r direction (dr)
 336   // in respect to angle phi within the ZR plane.
 337   // The histogramm has nx times ny entries.
 338   //
 339   TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
 340                      nz,-250.,250.,nr,85.,250.);
 341   Float_t x[3],dx[3];
 342   for (Int_t ir=1;ir<=nr;++ir) {
 343     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 344     x[0]=radius*TMath::Cos(phi);
 345     x[1]=radius*TMath::Sin(phi);
 346     for (Int_t iz=1;iz<=nz;++iz) {
 347       x[2]=h->GetXaxis()->GetBinCenter(iz);
 348       Int_t roc=x[2]>0.?0:18; // FIXME
 349       GetCorrection(x,roc,dx);
 350       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 351       Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 352       h->SetBinContent(iz,ir,r1-r0);
 353     }
 354   }
 355   return h;
 356
 357 }
 358
 359 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
 360   //
 361   // Simple plot functionality.
 362   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 363   // in respect to angle phi within the ZR plane.
 364   // The histogramm has nx times ny entries.
 365   //
 366   TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
 367                      nz,-250.,250.,nr,85.,250.);
 368   Float_t x[3],dx[3];
 369   for (Int_t iz=1;iz<=nz;++iz) {
 370     x[2]=h->GetXaxis()->GetBinCenter(iz);
 371     Int_t roc=x[2]>0.?0:18; // FIXME
 372     for (Int_t ir=1;ir<=nr;++ir) {
 373       Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 374       x[0]=radius*TMath::Cos(phi);
 375       x[1]=radius*TMath::Sin(phi);
 376       GetCorrection(x,roc,dx);
 377       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 378       Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 379       Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 380
 381       Float_t dphi=phi1-phi0;
 382       if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 383       if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 384
 385       h->SetBinContent(iz,ir,r0*dphi);
 386     }
 387   }
 388   return h;
 389 }
 390
 391 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
 392                                    const char *xlabel,const char *ylabel,const char *zlabel,
 393                                   Int_t nbinsx,Double_t xlow,Double_t xup,
 394                                   Int_t nbinsy,Double_t ylow,Double_t yup) {
 395   //
 396   // Helper function to create a 2d histogramm of given size
 397   //
 398
 399   TString hname=name;
 400   Int_t i=0;
 401   if (gDirectory) {
 402     while (gDirectory->FindObject(hname.Data())) {
 403       hname =name;
 404       hname+="_";
 405       hname+=i;
 406       ++i;
 407     }
 408   }
 409   TH2F *h=new TH2F(hname.Data(),title,
 410                    nbinsx,xlow,xup,
 411                    nbinsy,ylow,yup);
 412   h->GetXaxis()->SetTitle(xlabel);
 413   h->GetYaxis()->SetTitle(ylabel);
 414   h->GetZaxis()->SetTitle(zlabel);
 415   h->SetStats(0);
 416   return h;
 417 }
 418
 419
 420 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
 421
 422 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
 423                                                   const Double_t er[kNZ][kNR], Double_t &erValue ) {
 424   //
 425   // Interpolate table - 2D interpolation
 426   //
 427   Double_t saveEr[10] ;
 428
 429   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
 430   Search( kNR,   fgkRList,  r,   fKLow   ) ;
 431   if ( fJLow < 0 ) fJLow = 0 ;   // check if out of range
 432   if ( fKLow < 0 ) fKLow = 0 ;
 433   if ( fJLow + order  >=    kNZ - 1 ) fJLow =   kNZ - 1 - order ;
 434   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 435
 436   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
 437       saveEr[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
 438   }
 439   erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z )   ;
 440
 441 }
 442
 443
 444 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
 445                                        const Int_t order, const Double_t x ) {
 446   //
 447   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
 448   //
 449
 450   Double_t y ;
 451   if ( order == 2 ) {                // Quadratic Interpolation = 2
 452     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
 453     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
 454     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
 455   } else {                           // Linear Interpolation = 1
 456     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
 457   }
 458
 459   return (y);
 460
 461 }
 462
 463
 464 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
 465   //
 466   // Search an ordered table by starting at the most recently used point
 467   //
 468
 469   Long_t middle, high ;
 470   Int_t  ascend = 0, increment = 1 ;
 471
 472   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
 473
 474   if ( low < 0 || low > n-1 ) {
 475     low = -1 ; high = n ;
 476   } else {                                            // Ordered Search phase
 477     if ( (Int_t)( x >= xArray[low] ) == ascend )  {
 478       if ( low == n-1 ) return ;
 479       high = low + 1 ;
 480       while ( (Int_t)( x >= xArray[high] ) == ascend ) {
 481         low = high ;
 482         increment *= 2 ;
 483         high = low + increment ;
 484         if ( high > n-1 )  {  high = n ; break ;  }
 485       }
 486     } else {
 487       if ( low == 0 )  {  low = -1 ;  return ;  }
 488       high = low - 1 ;
 489       while ( (Int_t)( x < xArray[low] ) == ascend ) {
 490         high = low ;
 491         increment *= 2 ;
 492         if ( increment >= high )  {  low = -1 ;  break ;  }
 493         else  low = high - increment ;
 494       }
 495     }
 496   }
 497
 498   while ( (high-low) != 1 ) {                     // Binary Search Phase
 499     middle = ( high + low ) / 2 ;
 500     if ( (Int_t)( x >= xArray[middle] ) == ascend )
 501       low = middle ;
 502     else
 503       high = middle ;
 504   }
 505
 506   if ( x == xArray[n-1] ) low = n-2 ;
 507   if ( x == xArray[0]   ) low = 0 ;
 508
 509 }
 510
 511 void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, const TMatrixD &chargeDensity,
 512                                            TMatrixD &arrayErOverEz, const Int_t rows,
 513                                            const Int_t columns, const Int_t iterations ) {
 514   //
 515   // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
 516   //
 517   // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
 518   // boundary conditions on the first and last rows, and the first and last columns.  The remainder of the
 519   // array can be blank or contain a preliminary guess at the solution.  The Charge density matrix contains
 520   // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
 521   // you wish to solve Laplaces equation however it should not contain random numbers or you will get
 522   // random numbers back as a solution.
 523   // Poisson's equation is solved by iteratively relaxing the matrix to the final solution.  In order to
 524   // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
 525   // space.  First it solves the problem on a very sparse grid by skipping rows and columns in the original
 526   // matrix.  Then it doubles the number of points and solves the problem again.  Then it doubles the
 527   // number of points and solves the problem again.  This happens several times until the maximum number
 528   // of points has been included in the array.
 529   //
 530   // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
 531   // So rows == 2**M + 1 and columns == 2**N + 1.  The number of rows and columns can be different.
 532   //
 533   // Original code by Jim Thomas (STAR TPC Collaboration)
 534   //
 535
 536   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
 537
 538   const Float_t  gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
 539   const Float_t  gridSizeZ   =  fgkTPCZ0 / (columns-1) ;
 540   const Float_t  ratio       =  gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
 541
 542   TMatrixD  arrayEr(rows,columns) ;
 543   TMatrixD  arrayEz(rows,columns) ;
 544
 545   //Check that number of rows and columns is suitable for a binary expansion
 546
 547   if ( !IsPowerOfTwo(rows-1) ) {
 548     AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
 549     return;
 550   }
 551   if ( !IsPowerOfTwo(columns-1) ) {
 552     AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
 553     return;
 554   }
 555
 556   // Solve Poisson's equation in cylindrical coordinates by relaxation technique
 557   // Allow for different size grid spacing in R and Z directions
 558   // Use a binary expansion of the size of the matrix to speed up the solution of the problem
 559
 560   Int_t iOne = (rows-1)/4 ;
 561   Int_t jOne = (columns-1)/4 ;
 562   // Solve for N in 2**N, add one.
 563   Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
 564
 565   for ( Int_t count = 0 ; count < loops ; count++ ) {
 566     // Loop while the matrix expands & the resolution increases.
 567
 568     Float_t tempGridSizeR = gridSizeR * iOne ;
 569     Float_t tempRatio     = ratio * iOne * iOne / ( jOne * jOne ) ;
 570     Float_t tempFourth    = 1.0 / (2.0 + 2.0*tempRatio) ;
 571
 572     // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
 573     std::vector<float> coef1(rows) ;
 574     std::vector<float> coef2(rows) ;
 575
 576     for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
 577        Float_t radius = fgkIFCRadius + i*gridSizeR ;
 578       coef1[i] = 1.0 + tempGridSizeR/(2*radius);
 579       coef2[i] = 1.0 - tempGridSizeR/(2*radius);
 580     }
 581
 582     TMatrixD sumChargeDensity(rows,columns) ;
 583
 584     for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 585       Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
 586       for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 587         if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
 588         else {
 589           // Add up all enclosed charge density contributions within 1/2 unit in all directions
 590           Float_t weight = 0.0 ;
 591           Float_t sum    = 0.0 ;
 592           sumChargeDensity(i,j) = 0.0 ;
 593           for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
 594             for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
 595               if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
 596               else
 597                 weight = 1.0 ;
 598               // Note that this is cylindrical geometry
 599               sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
 600               sum += weight*radius ;
 601             }
 602           }
 603           sumChargeDensity(i,j) /= sum ;
 604         }
 605         sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
 606        }
 607     }
 608
 609     for ( Int_t k = 1 ; k <= iterations; k++ ) {
 610       // Solve Poisson's Equation
 611       // Over-relaxation index, must be >= 1 but < 2.  Arrange for it to evolve from 2 => 1
 612       // as interations increase.
 613       Float_t overRelax   = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
 614       Float_t overRelaxM1 = overRelax - 1.0 ;
 615       Float_t overRelaxtempFourth, overRelaxcoef5 ;
 616       overRelaxtempFourth = overRelax * tempFourth ;
 617       overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
 618
 619       for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 620         for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 621
 622           arrayV(i,j) = (   coef2[i]       *   arrayV(i-iOne,j)
 623                           + tempRatio      * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
 624                           - overRelaxcoef5 *   arrayV(i,j)
 625                           + coef1[i]       *   arrayV(i+iOne,j)
 626                           + sumChargeDensity(i,j)
 627                         ) * overRelaxtempFourth;
 628         }
 629       }
 630
 631       if ( k == iterations ) {
 632         // After full solution is achieved, copy low resolution solution into higher res array
 633         for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 634           for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 635
 636             if ( iOne > 1 ) {
 637               arrayV(i+iOne/2,j)                    =  ( arrayV(i+iOne,j) + arrayV(i,j)     ) / 2 ;
 638               if ( i == iOne )  arrayV(i-iOne/2,j) =  ( arrayV(0,j)       + arrayV(iOne,j) ) / 2 ;
 639             }
 640             if ( jOne > 1 ) {
 641               arrayV(i,j+jOne/2)                    =  ( arrayV(i,j+jOne) + arrayV(i,j) )     / 2 ;
 642               if ( j == jOne )  arrayV(i,j-jOne/2) =  ( arrayV(i,0)       + arrayV(i,jOne) ) / 2 ;
 643             }
 644             if ( iOne > 1 && jOne > 1 ) {
 645               arrayV(i+iOne/2,j+jOne/2) =  ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
 646               if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
 647               if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
 648               // Note that this leaves a point at the upper left and lower right corners uninitialized.
 649               // -> Not a big deal.
 650             }
 651
 652           }
 653         }
 654       }
 655
 656     }
 657
 658     iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
 659     jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
 660
 661   }
 662
 663   // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
 664   for ( Int_t j = 0 ; j < columns ; j++ ) {
 665     for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
 666     arrayEr(0,j)      =  -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
 667     arrayEr(rows-1,j) =  -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
 668   }
 669
 670   // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
 671   for ( Int_t i = 0 ; i < rows ; i++) {
 672     for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
 673     arrayEz(i,0)         =  -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
 674     arrayEz(i,columns-1) =  -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
 675   }
 676
 677   for ( Int_t i = 0 ; i < rows ; i++) {
 678     // Note: go back and compare to old version of this code.  See notes below.
 679     // JT Test ... attempt to divide by real Ez not Ez to first order
 680     for ( Int_t j = 0 ; j < columns ; j++ ) {
 681       arrayEz(i,j) += ezField;
 682       // This adds back the overall Z gradient of the field (main E field component)
 683     }
 684     // Warning: (-=) assumes you are using an error potetial without the overall Field included
 685   }
 686
 687   // Integrate Er/Ez from Z to zero
 688   for ( Int_t j = 0 ; j < columns ; j++ )  {
 689     for ( Int_t i = 0 ; i < rows ; i++ ) {
 690       Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
 691       arrayErOverEz(i,j) = 0.0 ;
 692       for ( Int_t k = j ; k < columns ; k++ ) {
 693         arrayErOverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
 694         if ( index != 4 )  index = 4; else index = 2 ;
 695       }
 696       if ( index == 4 ) arrayErOverEz(i,j)  -=  (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
 697       if ( index == 2 ) arrayErOverEz(i,j)  +=
 698         (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
 699                             -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) )   ;
 700       if ( j == columns-2 ) arrayErOverEz(i,j) =
 701         (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
 702                             +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
 703       if ( j == columns-1 ) arrayErOverEz(i,j) =  0.0 ;
 704     }
 705   }
 706
 707 }
 708
 709
 710
 711 Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
 712   //
 713   // Helperfunction: Check if integer is a power of 2
 714   //
 715   Int_t j = 0;
 716   while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
 717   if ( j == 1 ) return(1) ;  // True
 718   return(0) ;                // False
 719 }
 720
 721
 722 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
 723   //
 724   // Fit the track parameters - without and with distortion
 725   // 1. Space points in the TPC are simulated along the trajectory
 726   // 2. Space points distorted
 727   // 3. Fits  the non distorted and distroted track to the reference plane at refX
 728   // 4. For visualization and debugging  purposes the space points and tracks can be stored  in the tree - using the TTreeSRedirector functionality
 729   //
 730   // trackIn   - input track parameters
 731   // refX     - reference X to fit the track
 732   // dir      - direction - out=1 or in=-1
 733   // pcstream -  debug streamer to check the results
 734   //
 735   // see AliExternalTrackParam.h documentation:
 736   // track1.fP[0] - local y (rphi)
 737   // track1.fP[1] - z
 738   // track1.fP[2] - sinus of local inclination angle
 739   // track1.fP[3] - tangent of deep angle
 740   // track1.fP[4] - 1/pt
 741
 742   AliTPCROC * roc = AliTPCROC::Instance();
 743   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
 744   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
 745   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 746   const Double_t kMaxSnp = 0.85;
 747   const Double_t kSigmaY=0.1;
 748   const Double_t kSigmaZ=0.1;
 749   const Double_t kMaxR=500;
 750   const Double_t kMaxZ=500;
 751   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 752   Int_t npoints1=0;
 753   Int_t npoints2=0;
 754
 755   AliExternalTrackParam  track(trackIn); //
 756   // generate points
 757   AliTrackPointArray pointArray0(npoints0);
 758   AliTrackPointArray pointArray1(npoints0);
 759   Double_t xyz[3];
 760   if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
 761   //
 762   // simulate the track
 763   Int_t npoints=0;
 764   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
 765   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
 766     if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
 767     track.GetXYZ(xyz);
 768     xyz[0]+=gRandom->Gaus(0,0.00005);
 769     xyz[1]+=gRandom->Gaus(0,0.00005);
 770     xyz[2]+=gRandom->Gaus(0,0.00005);
 771     if (TMath::Abs(track.GetZ())>kMaxZ) break;
 772     if (TMath::Abs(track.GetX())>kMaxR) break;
 773     AliTrackPoint pIn0;                               // space point
 774     AliTrackPoint pIn1;
 775     Int_t sector= (xyz[2]>0)? 0:18;
 776     pointArray0.GetPoint(pIn0,npoints);
 777     pointArray1.GetPoint(pIn1,npoints);
 778     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 779     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
 780     DistortPoint(distPoint, sector);
 781     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
 782     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
 783     //
 784     track.Rotate(alpha);
 785     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
 786     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
 787     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
 788     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
 789     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
 790     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
 791     pointArray0.AddPoint(npoints, &pIn0);
 792     pointArray1.AddPoint(npoints, &pIn1);
 793     npoints++;
 794     if (npoints>=npoints0) break;
 795   }
 796   if (npoints<npoints0/2) return 0;
 797   //
 798   // refit track
 799   //
 800   AliExternalTrackParam *track0=0;
 801   AliExternalTrackParam *track1=0;
 802   AliTrackPoint   point1,point2,point3;
 803   if (dir==1) {  //make seed inner
 804     pointArray0.GetPoint(point1,1);
 805     pointArray0.GetPoint(point2,30);
 806     pointArray0.GetPoint(point3,60);
 807   }
 808   if (dir==-1){ //make seed outer
 809     pointArray0.GetPoint(point1,npoints-60);
 810     pointArray0.GetPoint(point2,npoints-30);
 811     pointArray0.GetPoint(point3,npoints-1);
 812   }
 813   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
 814   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
 815
 816   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
 817     Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
 818     //
 819     AliTrackPoint pIn0;
 820     AliTrackPoint pIn1;
 821     pointArray0.GetPoint(pIn0,ipoint);
 822     pointArray1.GetPoint(pIn1,ipoint);
 823     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
 824     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
 825     //
 826     if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
 827     if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
 828     if (TMath::Abs(track0->GetZ())>kMaxZ) break;
 829     if (TMath::Abs(track0->GetX())>kMaxR) break;
 830     if (TMath::Abs(track1->GetZ())>kMaxZ) break;
 831     if (TMath::Abs(track1->GetX())>kMaxR) break;
 832
 833     track.GetXYZ(xyz);  // distorted track also propagated to the same reference radius
 834     //
 835     Double_t pointPos[2]={0,0};
 836     Double_t pointCov[3]={0,0,0};
 837     pointPos[0]=prot0.GetY();//local y
 838     pointPos[1]=prot0.GetZ();//local z
 839     pointCov[0]=prot0.GetCov()[3];//simay^2
 840     pointCov[1]=prot0.GetCov()[4];//sigmayz
 841     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
 842     if (!track0->Update(pointPos,pointCov)) break;
 843     //
 844     Double_t deltaX=prot1.GetX()-prot0.GetX();   // delta X
 845     Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp()));  // deltaY due  delta X
 846     Double_t deltaZX=deltaX*track1->GetTgl();                           // deltaZ due  delta X
 847
 848     pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
 849     pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
 850     pointCov[0]=prot1.GetCov()[3];//simay^2
 851     pointCov[1]=prot1.GetCov()[4];//sigmayz
 852     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
 853     if (!track1->Update(pointPos,pointCov)) break;
 854     npoints1++;
 855     npoints2++;
 856   }
 857   if (npoints2<npoints)  return 0;
 858   AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
 859   track1->Rotate(track0->GetAlpha());
 860   AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
 861
 862   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
 863     "point0.="<<&pointArray0<<   //  points
 864     "point1.="<<&pointArray1<<   //  distorted points
 865     "trackIn.="<<&track<<       //  original track
 866     "track0.="<<track0<<         //  fitted track
 867     "track1.="<<track1<<         //  fitted distorted track
 868     "\n";
 869   new(&trackIn) AliExternalTrackParam(*track0);
 870   delete track0;
 871   return track1;
 872 }
 873
 874
 875
 876
 877
 878 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
 879   //
 880   // create the distortion tree on a mesh with granularity given by step
 881   // return the tree with distortions at given position
 882   // Map is created on the mesh with given step size
 883   //
 884   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
 885   Float_t xyz[3];
 886   for (Double_t x= -250; x<250; x+=step){
 887     for (Double_t y= -250; y<250; y+=step){
 888       Double_t r    = TMath::Sqrt(x*x+y*y);
 889       if (r<80) continue;
 890       if (r>250) continue;
 891       for (Double_t z= -250; z<250; z+=step){
 892         Int_t roc=(z>0)?0:18;
 893         xyz[0]=x;
 894         xyz[1]=y;
 895         xyz[2]=z;
 896         Double_t phi  = TMath::ATan2(y,x);
 897         DistortPoint(xyz,roc);
 898         Double_t r1    = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 899         Double_t phi1  = TMath::ATan2(xyz[1],xyz[0]);
 900         if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
 901         if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
 902         Double_t dx = xyz[0]-x;
 903         Double_t dy = xyz[1]-y;
 904         Double_t dz = xyz[2]-z;
 905         Double_t dr=r1-r;
 906         Double_t drphi=(phi1-phi)*r;
 907         (*pcstream)<<"distortion"<<
 908           "x="<<x<<           // original position
 909           "y="<<y<<
 910           "z="<<z<<
 911           "r="<<r<<
 912           "phi="<<phi<<
 913           "x1="<<xyz[0]<<      // distorted position
 914           "y1="<<xyz[1]<<
 915           "z1="<<xyz[2]<<
 916           "r1="<<r1<<
 917           "phi1="<<phi1<<
 918           //
 919           "dx="<<dx<<          // delta position
 920           "dy="<<dy<<
 921           "dz="<<dz<<
 922           "dr="<<dr<<
 923           "drphi="<<drphi<<
 924           "\n";
 925       }
 926     }
 927   }
 928   delete pcstream;
 929   TFile f(Form("correction%s.root",GetName()));
 930   TTree * tree = (TTree*)f.Get("distortion");
 931   TTree * tree2= tree->CopyTree("1");
 932   tree2->SetName(Form("dist%s",GetName()));
 933   tree2->SetDirectory(0);
 934   delete tree;
 935   return tree2;
 936 }
 937
 938
 939
 940
 941 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
 942   //
 943   // Make a fit tree:
 944   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
 945   // calculates partial distortions
 946   // Partial distortion is stored in the resulting tree
 947   // Output is storred in the file distortion_<dettype>_<partype>.root
 948   // Partial  distortion is stored with the name given by correction name
 949   //
 950   //
 951   // Parameters of function:
 952   // input     - input tree
 953   // dtype     - distortion type 0 - ITSTPC,  1 -TPCTRD, 2 - TPCvertex
 954   // ppype     - parameter type
 955   // corrArray - array with partial corrections
 956   // step      - skipe entries  - if 1 all entries processed - it is slow
 957   // debug     0 if debug on also space points dumped - it is slow
 958   const Double_t kMaxSnp = 0.85;
 959   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 960   //  const Double_t kB2C=-0.299792458e-3;
 961   const Int_t kMinEntries=50;
 962   Double_t phi,theta, snp, mean,rms, entries;
 963   tinput->SetBranchAddress("theta",&theta);
 964   tinput->SetBranchAddress("phi", &phi);
 965   tinput->SetBranchAddress("snp",&snp);
 966   tinput->SetBranchAddress("mean",&mean);
 967   tinput->SetBranchAddress("rms",&rms);
 968   tinput->SetBranchAddress("entries",&entries);
 969   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
 970   //
 971   Int_t nentries=tinput->GetEntries();
 972   Int_t ncorr=corrArray->GetEntries();
 973   Double_t corrections[100]={0}; //
 974   Double_t tPar[5];
 975   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 976   Double_t refX=0;
 977   Int_t dir=0;
 978   if (dtype==0) {refX=85.; dir=-1;}
 979   if (dtype==1) {refX=275.; dir=1;}
 980   if (dtype==2) {refX=85.; dir=-1;}
 981   if (dtype==3) {refX=360.; dir=-1;}
 982   //
 983   for (Int_t ientry=0; ientry<nentries; ientry+=step){
 984     tinput->GetEntry(ientry);
 985     if (TMath::Abs(snp)>kMaxSnp) continue;
 986     tPar[0]=0;
 987     tPar[1]=theta*refX;
 988     tPar[2]=snp;
 989     tPar[3]=theta;
 990     tPar[4]=(gRandom->Rndm()-0.5)*0.02;  // should be calculated - non equal to 0
 991     Double_t bz=AliTrackerBase::GetBz();
 992     if (refX>10. && TMath::Abs(bz)>0.1 )  tPar[4]=snp/(refX*bz*kB2C*2);
 993     tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
 994     AliExternalTrackParam track(refX,phi,tPar,cov);
 995     Double_t xyz[3];
 996     track.GetXYZ(xyz);
 997     Int_t id=0;
 998     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
 999     if (ptype==4 &&bz<0) mean*=-1;  // interpret as curvature
1000     (*pcstream)<<"fit"<<
1001       "bz="<<bz<<         // magnetic filed used
1002       "dtype="<<dtype<<   // detector match type
1003       "ptype="<<ptype<<   // parameter type
1004       "theta="<<theta<<   // theta
1005       "phi="<<phi<<       // phi
1006       "snp="<<snp<<       // snp
1007       "mean="<<mean<<     // mean dist value
1008       "rms="<<rms<<       // rms
1009       "gx="<<xyz[0]<<         // global position at reference
1010       "gy="<<xyz[1]<<         // global position at reference
1011       "gz="<<xyz[2]<<         // global position at reference
1012       "dRrec="<<dRrec<<      // delta Radius in reconstruction
1013       "id="<<id<<             // track id
1014       "entries="<<entries;// number of entries in bin
1015     //
1016     for (Int_t icorr=0; icorr<ncorr; icorr++) {
1017       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1018       corrections[icorr]=0;
1019       if (entries>kMinEntries){
1020         AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1021         AliExternalTrackParam *trackOut = 0;
1022         if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1023         if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1024         if (dtype==0) {refX=85.; dir=-1;}
1025         if (dtype==1) {refX=275.; dir=1;}
1026         if (dtype==2) {refX=0; dir=-1;}
1027         if (dtype==3) {refX=360.; dir=-1;}
1028         //
1029         if (trackOut){
1030           AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1031           trackOut->Rotate(trackIn.GetAlpha());
1032           trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1033           //
1034           corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1035           delete trackOut;
1036         }else{
1037           corrections[icorr]=0;
1038         }
1039         if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature
1040       }
1041       Double_t dRdummy=0;
1042       (*pcstream)<<"fit"<<
1043         Form("%s=",corr->GetName())<<corrections[icorr]<<   // dump correction value
1044         Form("dR%s=",corr->GetName())<<dRdummy;   // dump dummy correction value not needed for tracks
1045                                                   // for points it is neccessary
1046     }
1047     (*pcstream)<<"fit"<<"\n";
1048   }
1049   delete pcstream;
1050 }
1051
1052
1053
1054 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1055   //
1056   // Make a laser fit tree for global minimization
1057   //
1058   const Double_t cutErrY=0.1;
1059   const Double_t cutErrZ=0.1;
1060   const Double_t kEpsilon=0.00000001;
1061   TVectorD *vecdY=0;
1062   TVectorD *vecdZ=0;
1063   TVectorD *veceY=0;
1064   TVectorD *veceZ=0;
1065   AliTPCLaserTrack *ltr=0;
1066   AliTPCLaserTrack::LoadTracks();
1067   tree->SetBranchAddress("dY.",&vecdY);
1068   tree->SetBranchAddress("dZ.",&vecdZ);
1069   tree->SetBranchAddress("eY.",&veceY);
1070   tree->SetBranchAddress("eZ.",&veceZ);
1071   tree->SetBranchAddress("LTr.",&ltr);
1072   Int_t entries= tree->GetEntries();
1073   TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1074   Double_t bz=AliTrackerBase::GetBz();
1075   //
1076
1077   for (Int_t ientry=0; ientry<entries; ientry++){
1078     tree->GetEntry(ientry);
1079     if (!ltr->GetVecGX()){
1080       ltr->UpdatePoints();
1081     }
1082     TVectorD * delta= (itype==0)? vecdY:vecdZ;
1083     TVectorD * err= (itype==0)? veceY:veceZ;
1084
1085     for (Int_t irow=0; irow<159; irow++){
1086       Int_t nentries = 1000;
1087       if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1088       if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1089       Int_t dtype=4;
1090       Double_t phi   =(*ltr->GetVecPhi())[irow];
1091       Double_t theta =ltr->GetTgl();
1092       Double_t mean=delta->GetMatrixArray()[irow];
1093       Double_t gx=0,gy=0,gz=0;
1094       Double_t snp = (*ltr->GetVecP2())[irow];
1095       Double_t rms = 0.1+err->GetMatrixArray()[irow];
1096       gx = (*ltr->GetVecGX())[irow];
1097       gy = (*ltr->GetVecGY())[irow];
1098       gz = (*ltr->GetVecGZ())[irow];
1099       Int_t bundle= ltr->GetBundle();
1100       Double_t dRrec=0;
1101       //
1102       // get delta R used in reconstruction
1103       AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
1104       AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1105       const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1106       Double_t xyz0[3]={gx,gy,gz};
1107       Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1108       //
1109       // old ExB correction
1110       //
1111       if(recoParam&&recoParam->GetUseExBCorrection()) {
1112         Double_t xyz1[3]={gx,gy,gz};
1113         calib->GetExB()->Correct(xyz0,xyz1);
1114         Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1115         dRrec=oldR-newR;
1116       }
1117       if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1118         Float_t xyz1[3]={gx,gy,gz};
1119         Int_t sector=(gz>0)?0:18;
1120         correction->CorrectPoint(xyz1, sector);
1121         Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1122         dRrec=oldR-newR;
1123       }
1124
1125
1126       (*pcstream)<<"fit"<<
1127         "bz="<<bz<<         // magnetic filed used
1128         "dtype="<<dtype<<   // detector match type
1129         "ptype="<<itype<<   // parameter type
1130         "theta="<<theta<<   // theta
1131         "phi="<<phi<<       // phi
1132         "snp="<<snp<<       // snp
1133         "mean="<<mean<<     // mean dist value
1134         "rms="<<rms<<       // rms
1135         "gx="<<gx<<         // global position
1136         "gy="<<gy<<         // global position
1137         "gz="<<gz<<         // global position
1138         "dRrec="<<dRrec<<      // delta Radius in reconstruction
1139         "id="<<bundle<<     //bundle
1140         "entries="<<nentries;// number of entries in bin
1141       //
1142       //
1143       Double_t ky = TMath::Tan(TMath::ASin(snp));
1144       Int_t ncorr = corrArray->GetEntries();
1145       Double_t r0   = TMath::Sqrt(gx*gx+gy*gy);
1146       Double_t phi0 = TMath::ATan2(gy,gx);
1147       Double_t distortions[1000]={0};
1148       Double_t distortionsR[1000]={0};
1149       for (Int_t icorr=0; icorr<ncorr; icorr++) {
1150         AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1151         Float_t distPoint[3]={gx,gy,gz};
1152         Int_t sector= (gz>0)? 0:18;
1153         if (r0>80){
1154           corr->DistortPoint(distPoint, sector);
1155         }
1156         // Double_t value=distPoint[2]-gz;
1157         if (itype==0){
1158           Double_t r1   = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1159           Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1160           Double_t drphi= r0*(phi1-phi0);
1161           Double_t dr   = r1-r0;
1162           distortions[icorr]  = drphi-ky*dr;
1163           distortionsR[icorr] = dr;
1164         }
1165         (*pcstream)<<"fit"<<
1166           Form("%s=",corr->GetName())<<distortions[icorr]<<    // dump correction value
1167           Form("dR%s=",corr->GetName())<<distortionsR[icorr];   // dump correction R  value
1168       }
1169       (*pcstream)<<"fit"<<"\n";
1170     }
1171   }
1172   delete pcstream;
1173 }
1174
1175
1176
1177 void   AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
1178   //
1179   // make a distortion map out ou fthe residual histogram
1180   // Results are written to the debug streamer - pcstream
1181   // Parameters:
1182   //   his0       - input (4D) residual histogram
1183   //   pcstream   - file to write the tree
1184   //   run        - run number
1185   // marian.ivanov@cern.ch
1186   const Int_t kMinEntries=50;
1187   Int_t nbins1=his0->GetAxis(1)->GetNbins();
1188   Int_t first1=his0->GetAxis(1)->GetFirst();
1189   Int_t last1 =his0->GetAxis(1)->GetLast();
1190   //
1191   Double_t bz=AliTrackerBase::GetBz();
1192   Int_t idim[4]={0,1,2,3};
1193   for (Int_t ibin1=first1; ibin1<last1; ibin1++){   //axis 1 - theta
1194     //
1195     his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1196     Double_t       x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1197     THnSparse * his1 = his0->Projection(4,idim);  // projected histogram according range1
1198     Int_t nbins3     = his1->GetAxis(3)->GetNbins();
1199     Int_t first3     = his1->GetAxis(3)->GetFirst();
1200     Int_t last3      = his1->GetAxis(3)->GetLast();
1201     //
1202
1203     for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){   // axis 3 - local angle
1204       his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1205       Double_t      x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1206       if (ibin3<first3) {
1207         his1->GetAxis(3)->SetRangeUser(-1,1);
1208         x3=0;
1209       }
1210       THnSparse * his3= his1->Projection(4,idim);         //projected histogram according selection 3
1211       Int_t nbins2    = his3->GetAxis(2)->GetNbins();
1212       Int_t first2    = his3->GetAxis(2)->GetFirst();
1213       Int_t last2     = his3->GetAxis(2)->GetLast();
1214       //
1215       for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1216         his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1217         Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1218         TH1 * hisDelta = his3->Projection(0);
1219         //
1220         Double_t entries = hisDelta->GetEntries();
1221         Double_t mean=0, rms=0;
1222         if (entries>kMinEntries){
1223           mean    = hisDelta->GetMean();
1224           rms = hisDelta->GetRMS();
1225         }
1226         (*pcstream)<<hname<<
1227           "run="<<run<<
1228           "bz="<<bz<<
1229           "theta="<<x1<<
1230           "phi="<<x2<<
1231           "snp="<<x3<<
1232           "entries="<<entries<<
1233           "mean="<<mean<<
1234           "rms="<<rms<<
1235           "\n";
1236         delete hisDelta;
1237         printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1238       }
1239       delete his3;
1240     }
1241     delete his1;
1242   }
1243 }
1244
1245
1246
1247
1248
1249 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1250   //
1251   // Store object in the OCDB
1252   // By default the object is stored in the current directory
1253   // default comment consit of user name and the date
1254   //
1255   TString ocdbStorage="";
1256   ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1257   AliCDBMetaData *metaData= new AliCDBMetaData();
1258   metaData->SetObjectClassName("AliTPCCorrection");
1259   metaData->SetResponsible("Marian Ivanov");
1260   metaData->SetBeamPeriod(1);
1261   metaData->SetAliRootVersion("05-25-01"); //root version
1262   TString userName=gSystem->GetFromPipe("echo $USER");
1263   TString date=gSystem->GetFromPipe("date");
1264
1265   if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1266   if (comment) metaData->SetComment(comment);
1267   AliCDBId* id1=NULL;
1268   id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1269   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1270   gStorage->Put(this, (*id1), metaData);
1271 }
1272
1273
1274 void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
1275
1276   AliVertexerTracks *vertexer = new AliVertexerTracks(5);// 5kGaus
1277
1278   TObjArray   ATrk;               // Original Track array of Aside
1279   TObjArray   dATrk;              // Distorted Track array of A side
1280   UShort_t    *AId = new UShort_t[nTracks];      // A side Track ID
1281   TObjArray   CTrk;
1282   TObjArray   dCTrk;
1283   UShort_t    *CId = new UShort_t [nTracks];
1284   Int_t ID=0;
1285   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1286   TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
1287   fpt.SetParameters(7.24,0.120);
1288   fpt.SetNpx(10000);
1289   for(Int_t nt=0; nt<nTracks; nt++){
1290     Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
1291     Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
1292     Double_t pt = fpt.GetRandom();// momentum for f1
1293     Short_t sign=1;
1294     if(gRandom->Rndm() < 0.5){
1295       sign =1;
1296     }else{
1297       sign=-1;
1298     }
1299
1300     Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1301     Double_t pxyz[3];
1302     pxyz[0]=pt*TMath::Cos(phi);
1303     pxyz[1]=pt*TMath::Sin(phi);
1304     pxyz[2]=pt*TMath::Tan(theta);
1305     Double_t cv[21]={0};
1306     AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
1307
1308     Double_t refX=1.;
1309     Int_t dir=-1;
1310     AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir,  NULL);
1311     if (!td) continue;
1312     if (pcstream) (*pcstream)<<"track"<<
1313       "eta="<<eta<<
1314       "theta="<<theta<<
1315       "tOrig.="<<t<<
1316       "td.="<<td<<
1317       "\n";
1318     if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
1319       if (td){
1320         dATrk.AddLast(td);
1321         ATrk.AddLast(t);
1322         Int_t nn=ATrk.GetEntriesFast();
1323         AId[nn]=ID;
1324       }
1325     }else if(( eta<-0.07 )&&( eta>-etaCuts )){
1326       if (td){
1327         dCTrk.AddLast(td);
1328         CTrk.AddLast(t);
1329         Int_t nn=CTrk.GetEntriesFast();
1330         CId[nn]=ID;
1331       }
1332     }
1333     ID++;
1334   }// end of track loop
1335
1336   vertexer->SetTPCMode();
1337   vertexer->SetConstraintOff();
1338
1339   aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dATrk,AId));
1340   avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&ATrk,AId));
1341   cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dCTrk,CId));
1342   cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&CTrk,CId));
1343   if (pcstream) (*pcstream)<<"vertex"<<
1344     "x="<<orgVertex[0]<<
1345     "y="<<orgVertex[1]<<
1346     "z="<<orgVertex[2]<<
1347     "av.="<<&aV<<              // distorted vertex A side
1348     "cv.="<<&cV<<              // distroted vertex C side
1349     "avO.="<<&avOrg<<         // original vertex A side
1350     "cvO.="<<&cvOrg<<
1351     "\n";
1352   delete []AId;
1353   delete []CId;
1354 }
1355
1356 void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
1357   //
1358   // make correction available for visualization using
1359   // TFormula, TFX and TTree::Draw
1360   // important in order to check corrections and also compute dervied variables
1361   // e.g correction partial derivatives
1362   //
1363   // NOTE - class is not owner of correction
1364   //
1365   if (!fgVisualCorrection) fgVisualCorrection=new TObjArray;
1366   fgVisualCorrection->AddAt(corr, position);
1367 }
1368
1369
1370
1371 Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t localX, Double_t kZ, Int_t axisType, Int_t corrType){
1372   //
1373   // calculate the correction at given position - check the geffCorr
1374   //
1375   if (!fgVisualCorrection) return 0;
1376   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1377   if (!corr) return 0;
1378   Double_t phi=sector*TMath::Pi()/9.;
1379   Double_t gx = localX*TMath::Cos(phi);
1380   Double_t gy = localX*TMath::Sin(phi);
1381   Double_t gz = localX*kZ;
1382   Int_t nsector=(gz>0) ? 0:18;
1383   //
1384   //
1385   //
1386   Float_t distPoint[3]={gx,gy,gz};
1387   corr->DistortPoint(distPoint, nsector);
1388   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1389   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1390   Double_t phi0=TMath::ATan2(gy,gx);
1391   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1392   if (axisType==0) return r1-r0;
1393   if (axisType==1) return (phi1-phi0)*r0;
1394   if (axisType==2) return distPoint[2]-gz;
1395   return phi1-phi0;
1396 }
1397
1398 Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
1399   //
1400   // return correction at given x,y,z
1401   //
1402   if (!fgVisualCorrection) return 0;
1403   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1404   if (!corr) return 0;
1405   Double_t phi0= TMath::ATan2(gy,gx);
1406   Int_t nsector=(gz>0) ? 0:18;
1407   Float_t distPoint[3]={gx,gy,gz};
1408   corr->DistortPoint(distPoint, nsector);
1409   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1410   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1411   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1412   if (axisType==0) return r1-r0;
1413   if (axisType==1) return (phi1-phi0)*r0;
1414   if (axisType==2) return distPoint[2]-gz;
1415   return phi1-phi0;
1416 }