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 #include "AliTPCParamSR.h"
  55
  56 #include "AliTPCCorrection.h"
  57 #include "AliLog.h"
  58
  59 #include "AliExternalTrackParam.h"
  60 #include "AliTrackPointArray.h"
  61 #include "TDatabasePDG.h"
  62 #include "AliTrackerBase.h"
  63 #include "AliTPCROC.h"
  64 #include "THnSparse.h"
  65
  66 #include "AliTPCLaserTrack.h"
  67 #include "AliESDVertex.h"
  68 #include "AliVertexerTracks.h"
  69 #include "TDatabasePDG.h"
  70 #include "TF1.h"
  71 #include "TRandom.h"
  72
  73 #include "TDatabasePDG.h"
  74
  75 #include "AliTPCTransform.h"
  76 #include "AliTPCcalibDB.h"
  77 #include "AliTPCExB.h"
  78
  79 #include "AliTPCRecoParam.h"
  80
  81
  82 ClassImp(AliTPCCorrection)
  83
  84
  85 TObjArray *AliTPCCorrection::fgVisualCorrection=0;
  86 // instance of correction for visualization
  87
  88
  89 // FIXME: the following values should come from the database
  90 const Double_t AliTPCCorrection::fgkTPCZ0    = 249.7;     // nominal gating grid position
  91 const Double_t AliTPCCorrection::fgkIFCRadius=  83.5;     // radius which renders the "18 rod manifold" best -> compare calc. of Jim Thomas
  92 // compare gkIFCRadius=  83.05: Mean Radius of the Inner Field Cage ( 82.43 min,  83.70 max) (cm)
  93 const Double_t AliTPCCorrection::fgkOFCRadius= 254.5;     // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
  94 const Double_t AliTPCCorrection::fgkZOffSet  =   0.2;     // Offset from CE: calculate all distortions closer to CE as if at this point
  95 const Double_t AliTPCCorrection::fgkCathodeV = -100000.0; // Cathode Voltage (volts)
  96 const Double_t AliTPCCorrection::fgkGG       =     -70.0; // Gating Grid voltage (volts)
  97
  98 const Double_t  AliTPCCorrection::fgkdvdE = 0.0024; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
  99
 100 const Double_t AliTPCCorrection::fgkEM = -1.602176487e-19/9.10938215e-31; // charge/mass in [C/kg]
 101 const Double_t AliTPCCorrection::fgke0 = 8.854187817e-12;                 // vacuum permittivity [A·s/(V·m)]
 102
 103
 104 AliTPCCorrection::AliTPCCorrection()
 105   : TNamed("correction_unity","unity"),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
 106 {
 107   //
 108   // default constructor
 109   //
 110   if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
 111
 112   InitLookUpfulcrums();
 113
 114 }
 115
 116 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
 117 : TNamed(name,title),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
 118 {
 119   //
 120   // default constructor, that set the name and title
 121   //
 122   if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
 123
 124   InitLookUpfulcrums();
 125
 126 }
 127
 128 AliTPCCorrection::~AliTPCCorrection() {
 129   //
 130   // virtual destructor
 131   //
 132 }
 133
 134 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
 135   //
 136   // Corrects the initial coordinates x (cartesian coordinates)
 137   // according to the given effect (inherited classes)
 138   // roc represents the TPC read out chamber (offline numbering convention)
 139   //
 140   Float_t dx[3];
 141   GetCorrection(x,roc,dx);
 142   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 143 }
 144
 145 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 146   //
 147   // Corrects the initial coordinates x (cartesian coordinates) and stores the new
 148   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 149   // roc represents the TPC read out chamber (offline numbering convention)
 150   //
 151   Float_t dx[3];
 152   GetCorrection(x,roc,dx);
 153   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 154 }
 155
 156 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
 157   //
 158   // Distorts the initial coordinates x (cartesian coordinates)
 159   // according to the given effect (inherited classes)
 160   // roc represents the TPC read out chamber (offline numbering convention)
 161   //
 162   Float_t dx[3];
 163   GetDistortion(x,roc,dx);
 164   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 165 }
 166
 167 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 168   //
 169   // Distorts the initial coordinates x (cartesian coordinates) and stores the new
 170   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 171   // roc represents the TPC read out chamber (offline numbering convention)
 172   //
 173   Float_t dx[3];
 174   GetDistortion(x,roc,dx);
 175   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 176 }
 177
 178 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
 179   //
 180   // This function delivers the correction values dx in respect to the inital coordinates x
 181   // roc represents the TPC read out chamber (offline numbering convention)
 182   // Note: The dx is overwritten by the inherited effectice class ...
 183   //
 184   for (Int_t j=0;j<3;++j) { dx[j]=0.; }
 185 }
 186
 187 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
 188   //
 189   // This function delivers the distortion values dx in respect to the inital coordinates x
 190   // roc represents the TPC read out chamber (offline numbering convention)
 191   //
 192   GetCorrection(x,roc,dx);
 193   for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
 194 }
 195
 196 void AliTPCCorrection::Init() {
 197   //
 198   // Initialization funtion (not used at the moment)
 199   //
 200 }
 201
 202 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
 203   //
 204   // Update function
 205   //
 206 }
 207
 208 void AliTPCCorrection::Print(Option_t* /*option*/) const {
 209   //
 210   // Print function to check which correction classes are used
 211   // option=="d" prints details regarding the setted magnitude
 212   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 213   //
 214   printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
 215 }
 216
 217 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
 218   //
 219   // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
 220   // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
 221   // calibration run
 222   //
 223   fT1=t1;
 224   fT2=t2;
 225   //SetOmegaTauT1T2(omegaTau, t1, t2);
 226 }
 227
 228 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
 229   //
 230   // Simple plot functionality.
 231   // Returns a 2d hisogram which represents the corrections in radial direction (dr)
 232   // in respect to position z within the XY plane.
 233   // The histogramm has nx times ny entries.
 234   //
 235   AliTPCParam* tpcparam = new AliTPCParamSR;
 236
 237   TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
 238                      nx,-250.,250.,ny,-250.,250.);
 239   Float_t x[3],dx[3];
 240   x[2]=z;
 241   Int_t roc=z>0.?0:18; // FIXME
 242   for (Int_t iy=1;iy<=ny;++iy) {
 243     x[1]=h->GetYaxis()->GetBinCenter(iy);
 244     for (Int_t ix=1;ix<=nx;++ix) {
 245       x[0]=h->GetXaxis()->GetBinCenter(ix);
 246       GetCorrection(x,roc,dx);
 247       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 248       if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
 249         Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 250         h->SetBinContent(ix,iy,r1-r0);
 251       }
 252       else
 253         h->SetBinContent(ix,iy,0.);
 254     }
 255   }
 256   delete tpcparam;
 257   return h;
 258 }
 259
 260 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
 261   //
 262   // Simple plot functionality.
 263   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 264   // in respect to position z within the XY plane.
 265   // The histogramm has nx times ny entries.
 266   //
 267
 268   AliTPCParam* tpcparam = new AliTPCParamSR;
 269
 270   TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
 271                      nx,-250.,250.,ny,-250.,250.);
 272   Float_t x[3],dx[3];
 273   x[2]=z;
 274   Int_t roc=z>0.?0:18; // FIXME
 275   for (Int_t iy=1;iy<=ny;++iy) {
 276     x[1]=h->GetYaxis()->GetBinCenter(iy);
 277     for (Int_t ix=1;ix<=nx;++ix) {
 278       x[0]=h->GetXaxis()->GetBinCenter(ix);
 279       GetCorrection(x,roc,dx);
 280       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 281       if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
 282         Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 283         Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 284
 285         Float_t dphi=phi1-phi0;
 286         if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 287         if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 288
 289         h->SetBinContent(ix,iy,r0*dphi);
 290       }
 291       else
 292         h->SetBinContent(ix,iy,0.);
 293     }
 294   }
 295   delete tpcparam;
 296   return h;
 297 }
 298
 299 TH2F* AliTPCCorrection::CreateHistoDZinXY(Float_t z,Int_t nx,Int_t ny) {
 300   //
 301   // Simple plot functionality.
 302   // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
 303   // in respect to position z within the XY plane.
 304   // The histogramm has nx times ny entries.
 305   //
 306
 307   AliTPCParam* tpcparam = new AliTPCParamSR;
 308
 309   TH2F *h=CreateTH2F("dz_xy",GetTitle(),"x [cm]","y [cm]","dz [cm]",
 310                      nx,-250.,250.,ny,-250.,250.);
 311   Float_t x[3],dx[3];
 312   x[2]=z;
 313   Int_t roc=z>0.?0:18; // FIXME
 314   for (Int_t iy=1;iy<=ny;++iy) {
 315     x[1]=h->GetYaxis()->GetBinCenter(iy);
 316     for (Int_t ix=1;ix<=nx;++ix) {
 317       x[0]=h->GetXaxis()->GetBinCenter(ix);
 318       GetCorrection(x,roc,dx);
 319       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 320       if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
 321         h->SetBinContent(ix,iy,dx[2]);
 322       }
 323       else
 324         h->SetBinContent(ix,iy,0.);
 325     }
 326   }
 327   delete tpcparam;
 328   return h;
 329 }
 330
 331 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
 332   //
 333   // Simple plot functionality.
 334   // Returns a 2d hisogram which represents the corrections in r direction (dr)
 335   // in respect to angle phi within the ZR plane.
 336   // The histogramm has nx times ny entries.
 337   //
 338   TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
 339                      nz,-250.,250.,nr,85.,250.);
 340   Float_t x[3],dx[3];
 341   for (Int_t ir=1;ir<=nr;++ir) {
 342     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 343     x[0]=radius*TMath::Cos(phi);
 344     x[1]=radius*TMath::Sin(phi);
 345     for (Int_t iz=1;iz<=nz;++iz) {
 346       x[2]=h->GetXaxis()->GetBinCenter(iz);
 347       Int_t roc=x[2]>0.?0:18; // FIXME
 348       GetCorrection(x,roc,dx);
 349       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 350       Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 351       h->SetBinContent(iz,ir,r1-r0);
 352     }
 353   }
 354   return h;
 355
 356 }
 357
 358 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
 359   //
 360   // Simple plot functionality.
 361   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 362   // in respect to angle phi within the ZR plane.
 363   // The histogramm has nx times ny entries.
 364   //
 365   TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
 366                      nz,-250.,250.,nr,85.,250.);
 367   Float_t x[3],dx[3];
 368   for (Int_t iz=1;iz<=nz;++iz) {
 369     x[2]=h->GetXaxis()->GetBinCenter(iz);
 370     Int_t roc=x[2]>0.?0:18; // FIXME
 371     for (Int_t ir=1;ir<=nr;++ir) {
 372       Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 373       x[0]=radius*TMath::Cos(phi);
 374       x[1]=radius*TMath::Sin(phi);
 375       GetCorrection(x,roc,dx);
 376       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 377       Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 378       Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 379
 380       Float_t dphi=phi1-phi0;
 381       if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 382       if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 383
 384       h->SetBinContent(iz,ir,r0*dphi);
 385     }
 386   }
 387   return h;
 388 }
 389
 390 TH2F* AliTPCCorrection::CreateHistoDZinZR(Float_t phi,Int_t nz,Int_t nr) {
 391   //
 392   // Simple plot functionality.
 393   // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
 394   // in respect to angle phi within the ZR plane.
 395   // The histogramm has nx times ny entries.
 396   //
 397   TH2F *h=CreateTH2F("dz_zr",GetTitle(),"z [cm]","r [cm]","dz [cm]",
 398                      nz,-250.,250.,nr,85.,250.);
 399   Float_t x[3],dx[3];
 400   for (Int_t ir=1;ir<=nr;++ir) {
 401     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 402     x[0]=radius*TMath::Cos(phi);
 403     x[1]=radius*TMath::Sin(phi);
 404     for (Int_t iz=1;iz<=nz;++iz) {
 405       x[2]=h->GetXaxis()->GetBinCenter(iz);
 406       Int_t roc=x[2]>0.?0:18; // FIXME
 407       GetCorrection(x,roc,dx);
 408       h->SetBinContent(iz,ir,dx[2]);
 409     }
 410   }
 411   return h;
 412
 413 }
 414
 415
 416 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
 417                                    const char *xlabel,const char *ylabel,const char *zlabel,
 418                                   Int_t nbinsx,Double_t xlow,Double_t xup,
 419                                   Int_t nbinsy,Double_t ylow,Double_t yup) {
 420   //
 421   // Helper function to create a 2d histogramm of given size
 422   //
 423
 424   TString hname=name;
 425   Int_t i=0;
 426   if (gDirectory) {
 427     while (gDirectory->FindObject(hname.Data())) {
 428       hname =name;
 429       hname+="_";
 430       hname+=i;
 431       ++i;
 432     }
 433   }
 434   TH2F *h=new TH2F(hname.Data(),title,
 435                    nbinsx,xlow,xup,
 436                    nbinsy,ylow,yup);
 437   h->GetXaxis()->SetTitle(xlabel);
 438   h->GetYaxis()->SetTitle(ylabel);
 439   h->GetZaxis()->SetTitle(zlabel);
 440   h->SetStats(0);
 441   return h;
 442 }
 443
 444 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
 445
 446 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
 447                                                   const Double_t er[kNZ][kNR], Double_t &erValue ) {
 448   //
 449   // Interpolate table - 2D interpolation
 450   //
 451   Double_t saveEr[5] = {0,0,0,0,0};
 452
 453   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
 454   Search( kNR,   fgkRList,  r,   fKLow   ) ;
 455   if ( fJLow < 0 ) fJLow = 0 ;   // check if out of range
 456   if ( fKLow < 0 ) fKLow = 0 ;
 457   if ( fJLow + order  >=    kNZ - 1 ) fJLow =   kNZ - 1 - order ;
 458   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 459
 460   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
 461       saveEr[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
 462   }
 463   erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z )   ;
 464
 465 }
 466
 467 void AliTPCCorrection::Interpolate3DEdistortion( const Int_t order, const Double_t r, const Float_t phi, const Double_t z,
 468                                                  const Double_t er[kNZ][kNPhi][kNR], const Double_t ephi[kNZ][kNPhi][kNR], const Double_t ez[kNZ][kNPhi][kNR],
 469                                                  Double_t &erValue, Double_t &ephiValue, Double_t &ezValue) {
 470   //
 471   // Interpolate table - 3D interpolation
 472   //
 473
 474   Double_t saveEr[5]= {0,0,0,0,0};
 475   Double_t savedEr[5]= {0,0,0,0,0} ;
 476
 477   Double_t saveEphi[5]= {0,0,0,0,0};
 478   Double_t savedEphi[5]= {0,0,0,0,0} ;
 479
 480   Double_t saveEz[5]= {0,0,0,0,0};
 481   Double_t savedEz[5]= {0,0,0,0,0} ;
 482
 483   Search( kNZ,   fgkZList,   z,   fILow   ) ;
 484   Search( kNPhi, fgkPhiList, z,   fJLow   ) ;
 485   Search( kNR,   fgkRList,   r,   fKLow   ) ;
 486
 487   if ( fILow < 0 ) fILow = 0 ;   // check if out of range
 488   if ( fJLow < 0 ) fJLow = 0 ;
 489   if ( fKLow < 0 ) fKLow = 0 ;
 490
 491   if ( fILow + order  >=    kNZ - 1 ) fILow =   kNZ - 1 - order ;
 492   if ( fJLow + order  >=  kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
 493   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 494
 495   for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
 496     for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
 497       saveEr[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r )   ;
 498       saveEphi[j-fJLow]   = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
 499       saveEz[j-fJLow]     = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r )   ;
 500     }
 501     savedEr[i-fILow]     = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi )   ;
 502     savedEphi[i-fILow]   = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
 503     savedEz[i-fILow]     = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi )   ;
 504   }
 505   erValue     = Interpolate( &fgkZList[fILow], savedEr, order, z )    ;
 506   ephiValue   = Interpolate( &fgkZList[fILow], savedEphi, order, z )  ;
 507   ezValue     = Interpolate( &fgkZList[fILow], savedEz, order, z )    ;
 508
 509 }
 510
 511 Double_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
 512                                               const Int_t nx,  const Int_t ny, const Double_t xv[], const Double_t yv[],
 513                                               const TMatrixD &array ) {
 514   //
 515   // Interpolate table (TMatrix format) - 2D interpolation
 516   //
 517
 518   static  Int_t jlow = 0, klow = 0 ;
 519   Double_t saveArray[5] = {0,0,0,0,0} ;
 520
 521   Search( nx,  xv,  x,   jlow  ) ;
 522   Search( ny,  yv,  y,   klow  ) ;
 523   if ( jlow < 0 ) jlow = 0 ;   // check if out of range
 524   if ( klow < 0 ) klow = 0 ;
 525   if ( jlow + order  >=    nx - 1 ) jlow =   nx - 1 - order ;
 526   if ( klow + order  >=    ny - 1 ) klow =   ny - 1 - order ;
 527
 528   for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
 529     {
 530       Double_t *ajkl = &((TMatrixD&)array)(j,klow);
 531       saveArray[j-jlow]  = Interpolate( &yv[klow], ajkl , order, y )   ;
 532     }
 533
 534   return( Interpolate( &xv[jlow], saveArray, order, x ) )   ;
 535
 536 }
 537
 538 Double_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x,   const Double_t y,   const Double_t z,
 539                                               const Int_t  nx,    const Int_t  ny,    const Int_t  nz,
 540                                               const Double_t xv[], const Double_t yv[], const Double_t zv[],
 541                                               TMatrixD **arrayofArrays ) {
 542   //
 543   // Interpolate table (TMatrix format) - 3D interpolation
 544   //
 545
 546   static  Int_t ilow = 0, jlow = 0, klow = 0 ;
 547   Double_t saveArray[5]= {0,0,0,0,0};
 548   Double_t savedArray[5]= {0,0,0,0,0} ;
 549
 550   Search( nx, xv, x, ilow   ) ;
 551   Search( ny, yv, y, jlow   ) ;
 552   Search( nz, zv, z, klow   ) ;
 553
 554   if ( ilow < 0 ) ilow = 0 ;   // check if out of range
 555   if ( jlow < 0 ) jlow = 0 ;
 556   if ( klow < 0 ) klow = 0 ;
 557
 558   if ( ilow + order  >=    nx - 1 ) ilow =   nx - 1 - order ;
 559   if ( jlow + order  >=    ny - 1 ) jlow =   ny - 1 - order ;
 560   if ( klow + order  >=    nz - 1 ) klow =   nz - 1 - order ;
 561
 562   for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
 563     {
 564       TMatrixD &table = *arrayofArrays[k] ;
 565       for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
 566         {
 567           saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y )   ;
 568         }
 569       savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x )   ;
 570     }
 571   return( Interpolate( &zv[klow], savedArray, order, z ) )   ;
 572
 573 }
 574
 575
 576 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
 577                                        const Int_t order, const Double_t x ) {
 578   //
 579   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
 580   //
 581
 582   Double_t y ;
 583   if ( order == 2 ) {                // Quadratic Interpolation = 2
 584     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
 585     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
 586     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
 587   } else {                           // Linear Interpolation = 1
 588     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
 589   }
 590
 591   return (y);
 592
 593 }
 594
 595
 596 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
 597   //
 598   // Search an ordered table by starting at the most recently used point
 599   //
 600
 601   Long_t middle, high ;
 602   Int_t  ascend = 0, increment = 1 ;
 603
 604   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
 605
 606   if ( low < 0 || low > n-1 ) {
 607     low = -1 ; high = n ;
 608   } else {                                            // Ordered Search phase
 609     if ( (Int_t)( x >= xArray[low] ) == ascend )  {
 610       if ( low == n-1 ) return ;
 611       high = low + 1 ;
 612       while ( (Int_t)( x >= xArray[high] ) == ascend ) {
 613         low = high ;
 614         increment *= 2 ;
 615         high = low + increment ;
 616         if ( high > n-1 )  {  high = n ; break ;  }
 617       }
 618     } else {
 619       if ( low == 0 )  {  low = -1 ;  return ;  }
 620       high = low - 1 ;
 621       while ( (Int_t)( x < xArray[low] ) == ascend ) {
 622         high = low ;
 623         increment *= 2 ;
 624         if ( increment >= high )  {  low = -1 ;  break ;  }
 625         else  low = high - increment ;
 626       }
 627     }
 628   }
 629
 630   while ( (high-low) != 1 ) {                     // Binary Search Phase
 631     middle = ( high + low ) / 2 ;
 632     if ( (Int_t)( x >= xArray[middle] ) == ascend )
 633       low = middle ;
 634     else
 635       high = middle ;
 636   }
 637
 638   if ( x == xArray[n-1] ) low = n-2 ;
 639   if ( x == xArray[0]   ) low = 0 ;
 640
 641 }
 642
 643 void AliTPCCorrection::InitLookUpfulcrums() {
 644   //
 645   // Initialization of interpolation points - for main look up table
 646   //   (course grid in the middle, fine grid on the borders)
 647   //
 648
 649   AliTPCROC * roc = AliTPCROC::Instance();
 650   const Double_t rLow =  TMath::Floor(roc->GetPadRowRadii(0,0))-1; // first padRow plus some margin
 651
 652   // fulcrums in R
 653   fgkRList[0] = rLow;
 654   for (Int_t i = 1; i<kNR; i++) {
 655     fgkRList[i] = fgkRList[i-1] + 3.5;     // 3.5 cm spacing
 656     if (fgkRList[i]<90 ||fgkRList[i]>245)
 657        fgkRList[i] = fgkRList[i-1] + 0.5; // 0.5 cm spacing
 658     else if (fgkRList[i]<100 || fgkRList[i]>235)
 659        fgkRList[i] = fgkRList[i-1] + 1.5;  // 1.5 cm spacing
 660     else if (fgkRList[i]<120 || fgkRList[i]>225)
 661        fgkRList[i] = fgkRList[i-1] + 2.5;  // 2.5 cm spacing
 662   }
 663
 664   // fulcrums in Z
 665   fgkZList[0] = -249.5;
 666   fgkZList[kNZ-1] = 249.5;
 667   for (Int_t j = 1; j<kNZ/2; j++) {
 668     fgkZList[j] = fgkZList[j-1];
 669     if      (TMath::Abs(fgkZList[j])< 0.15)
 670       fgkZList[j] = fgkZList[j-1] + 0.09; // 0.09 cm spacing
 671     else if(TMath::Abs(fgkZList[j])< 0.6)
 672       fgkZList[j] = fgkZList[j-1] + 0.4; // 0.4 cm spacing
 673     else if      (TMath::Abs(fgkZList[j])< 2.5 || TMath::Abs(fgkZList[j])>248)
 674       fgkZList[j] = fgkZList[j-1] + 0.5; // 0.5 cm spacing
 675     else if (TMath::Abs(fgkZList[j])<10 || TMath::Abs(fgkZList[j])>235)
 676       fgkZList[j] = fgkZList[j-1] + 1.5;  // 1.5 cm spacing
 677     else if (TMath::Abs(fgkZList[j])<25 || TMath::Abs(fgkZList[j])>225)
 678       fgkZList[j] = fgkZList[j-1] + 2.5;  // 2.5 cm spacing
 679     else
 680       fgkZList[j] = fgkZList[j-1] + 4;  // 4 cm spacing
 681
 682     fgkZList[kNZ-j-1] = -fgkZList[j];
 683   }
 684
 685   // fulcrums in phi
 686   for (Int_t k = 0; k<kNPhi; k++)
 687     fgkPhiList[k] = TMath::TwoPi()*k/(kNPhi-1);
 688
 689
 690 }
 691
 692
 693 void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
 694                                            TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
 695                                            const Int_t rows, const Int_t columns, const Int_t iterations,
 696                                            const Bool_t rocDisplacement ) {
 697   //
 698   // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
 699   //
 700   // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
 701   // boundary conditions on the first and last rows, and the first and last columns.  The remainder of the
 702   // array can be blank or contain a preliminary guess at the solution.  The Charge density matrix contains
 703   // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
 704   // you wish to solve Laplaces equation however it should not contain random numbers or you will get
 705   // random numbers back as a solution.
 706   // Poisson's equation is solved by iteratively relaxing the matrix to the final solution.  In order to
 707   // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
 708   // space.  First it solves the problem on a very sparse grid by skipping rows and columns in the original
 709   // matrix.  Then it doubles the number of points and solves the problem again.  Then it doubles the
 710   // number of points and solves the problem again.  This happens several times until the maximum number
 711   // of points has been included in the array.
 712   //
 713   // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
 714   // So rows == 2**M + 1 and columns == 2**N + 1.  The number of rows and columns can be different.
 715   //
 716   // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
 717   //
 718   // Original code by Jim Thomas (STAR TPC Collaboration)
 719   //
 720
 721   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
 722
 723   const Float_t  gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
 724   const Float_t  gridSizeZ   =  fgkTPCZ0 / (columns-1) ;
 725   const Float_t  ratio       =  gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
 726
 727   TMatrixD  arrayEr(rows,columns) ;
 728   TMatrixD  arrayEz(rows,columns) ;
 729
 730   //Check that number of rows and columns is suitable for a binary expansion
 731
 732   if ( !IsPowerOfTwo(rows-1) ) {
 733     AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
 734     return;
 735   }
 736   if ( !IsPowerOfTwo(columns-1) ) {
 737     AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
 738     return;
 739   }
 740
 741   // Solve Poisson's equation in cylindrical coordinates by relaxation technique
 742   // Allow for different size grid spacing in R and Z directions
 743   // Use a binary expansion of the size of the matrix to speed up the solution of the problem
 744
 745   Int_t iOne = (rows-1)/4 ;
 746   Int_t jOne = (columns-1)/4 ;
 747   // Solve for N in 2**N, add one.
 748   Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
 749
 750   for ( Int_t count = 0 ; count < loops ; count++ ) {
 751     // Loop while the matrix expands & the resolution increases.
 752
 753     Float_t tempGridSizeR = gridSizeR * iOne ;
 754     Float_t tempRatio     = ratio * iOne * iOne / ( jOne * jOne ) ;
 755     Float_t tempFourth    = 1.0 / (2.0 + 2.0*tempRatio) ;
 756
 757     // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
 758     std::vector<float> coef1(rows) ;
 759     std::vector<float> coef2(rows) ;
 760
 761     for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
 762        Float_t radius = fgkIFCRadius + i*gridSizeR ;
 763       coef1[i] = 1.0 + tempGridSizeR/(2*radius);
 764       coef2[i] = 1.0 - tempGridSizeR/(2*radius);
 765     }
 766
 767     TMatrixD sumChargeDensity(rows,columns) ;
 768
 769     for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 770       Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
 771       for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 772         if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
 773         else {
 774           // Add up all enclosed charge density contributions within 1/2 unit in all directions
 775           Float_t weight = 0.0 ;
 776           Float_t sum    = 0.0 ;
 777           sumChargeDensity(i,j) = 0.0 ;
 778           for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
 779             for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
 780               if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
 781               else
 782                 weight = 1.0 ;
 783               // Note that this is cylindrical geometry
 784               sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
 785               sum += weight*radius ;
 786             }
 787           }
 788           sumChargeDensity(i,j) /= sum ;
 789         }
 790         sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
 791        }
 792     }
 793
 794     for ( Int_t k = 1 ; k <= iterations; k++ ) {
 795       // Solve Poisson's Equation
 796       // Over-relaxation index, must be >= 1 but < 2.  Arrange for it to evolve from 2 => 1
 797       // as interations increase.
 798       Float_t overRelax   = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
 799       Float_t overRelaxM1 = overRelax - 1.0 ;
 800       Float_t overRelaxtempFourth, overRelaxcoef5 ;
 801       overRelaxtempFourth = overRelax * tempFourth ;
 802       overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
 803
 804       for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 805         for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 806
 807           arrayV(i,j) = (   coef2[i]       *   arrayV(i-iOne,j)
 808                           + tempRatio      * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
 809                           - overRelaxcoef5 *   arrayV(i,j)
 810                           + coef1[i]       *   arrayV(i+iOne,j)
 811                           + sumChargeDensity(i,j)
 812                         ) * overRelaxtempFourth;
 813         }
 814       }
 815
 816       if ( k == iterations ) {
 817         // After full solution is achieved, copy low resolution solution into higher res array
 818         for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 819           for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 820
 821             if ( iOne > 1 ) {
 822               arrayV(i+iOne/2,j)                    =  ( arrayV(i+iOne,j) + arrayV(i,j)     ) / 2 ;
 823               if ( i == iOne )  arrayV(i-iOne/2,j) =  ( arrayV(0,j)       + arrayV(iOne,j) ) / 2 ;
 824             }
 825             if ( jOne > 1 ) {
 826               arrayV(i,j+jOne/2)                    =  ( arrayV(i,j+jOne) + arrayV(i,j) )     / 2 ;
 827               if ( j == jOne )  arrayV(i,j-jOne/2) =  ( arrayV(i,0)       + arrayV(i,jOne) ) / 2 ;
 828             }
 829             if ( iOne > 1 && jOne > 1 ) {
 830               arrayV(i+iOne/2,j+jOne/2) =  ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
 831               if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
 832               if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
 833               // Note that this leaves a point at the upper left and lower right corners uninitialized.
 834               // -> Not a big deal.
 835             }
 836
 837           }
 838         }
 839       }
 840
 841     }
 842
 843     iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
 844     jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
 845
 846     sumChargeDensity.Clear();
 847   }
 848
 849   // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
 850   for ( Int_t j = 0 ; j < columns ; j++ ) {
 851     for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
 852     arrayEr(0,j)      =  -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
 853     arrayEr(rows-1,j) =  -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
 854   }
 855
 856   // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
 857   for ( Int_t i = 0 ; i < rows ; i++) {
 858     for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
 859     arrayEz(i,0)         =  -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
 860     arrayEz(i,columns-1) =  -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
 861   }
 862
 863   for ( Int_t i = 0 ; i < rows ; i++) {
 864     // Note: go back and compare to old version of this code.  See notes below.
 865     // JT Test ... attempt to divide by real Ez not Ez to first order
 866     for ( Int_t j = 0 ; j < columns ; j++ ) {
 867       arrayEz(i,j) += ezField;
 868       // This adds back the overall Z gradient of the field (main E field component)
 869     }
 870     // Warning: (-=) assumes you are using an error potetial without the overall Field included
 871   }
 872
 873   // Integrate Er/Ez from Z to zero
 874   for ( Int_t j = 0 ; j < columns ; j++ )  {
 875     for ( Int_t i = 0 ; i < rows ; i++ ) {
 876
 877       Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
 878       arrayErOverEz(i,j) = 0.0 ;
 879       arrayDeltaEz(i,j) = 0.0 ;
 880
 881       for ( Int_t k = j ; k < columns ; k++ ) {
 882         arrayErOverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
 883         arrayDeltaEz(i,j)   +=  index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
 884         if ( index != 4 )  index = 4; else index = 2 ;
 885       }
 886       if ( index == 4 ) {
 887         arrayErOverEz(i,j)  -=  (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
 888         arrayDeltaEz(i,j)   -=  (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
 889       }
 890       if ( index == 2 ) {
 891         arrayErOverEz(i,j)  +=  (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
 892                                                     -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
 893         arrayDeltaEz(i,j)   +=  (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
 894                                                     -2.5*(arrayEz(i,columns-1)-ezField));
 895       }
 896       if ( j == columns-2 ) {
 897         arrayErOverEz(i,j) =  (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
 898                                                   +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
 899         arrayDeltaEz(i,j)  =  (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
 900                                                   +1.5*(arrayEz(i,columns-1)-ezField) ) ;
 901       }
 902       if ( j == columns-1 ) {
 903         arrayErOverEz(i,j) =  0.0 ;
 904         arrayDeltaEz(i,j)  =  0.0 ;
 905       }
 906     }
 907   }
 908
 909   // calculate z distortion from the integrated Delta Ez residuals
 910   // and include the aquivalence (Volt to cm) of the ROC shift !!
 911
 912   for ( Int_t j = 0 ; j < columns ; j++ )  {
 913     for ( Int_t i = 0 ; i < rows ; i++ ) {
 914
 915       // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
 916       arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
 917
 918       // ROC Potential in cm aquivalent
 919       Double_t dzROCShift =  arrayV(i, columns -1)/ezField;
 920       if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift;  // add the ROC misaligment
 921
 922     }
 923   }
 924
 925   arrayEr.Clear();
 926   arrayEz.Clear();
 927
 928 }
 929
 930 void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
 931                     TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
 932                     const Int_t rows, const Int_t columns,  const Int_t phislices,
 933                     const Float_t deltaphi, const Int_t iterations, const Int_t symmetry,
 934                     Bool_t rocDisplacement  ) {
 935   //
 936   // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
 937   //
 938   //    NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
 939   //    The number of rows and COLUMNS can be different.
 940   //
 941   //    ROWS       ==  2**M + 1
 942   //    COLUMNS    ==  2**N + 1
 943   //    PHISLICES  ==  Arbitrary but greater than 3
 944   //
 945   //    DeltaPhi in Radians
 946   //
 947   //    SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
 948   //             = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
 949   //
 950   // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
 951
 952   const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
 953
 954   const Float_t  gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
 955   const Float_t  gridSizePhi =  deltaphi ;
 956   const Float_t  gridSizeZ   =  fgkTPCZ0 / (columns-1) ;
 957   const Float_t  ratioPhi    =  gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
 958   const Float_t  ratioZ      =  gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
 959
 960   TMatrixD arrayE(rows,columns) ;
 961
 962   // Check that the number of rows and columns is suitable for a binary expansion
 963   if ( !IsPowerOfTwo((rows-1))    ) {
 964     AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
 965     return; }
 966   if ( !IsPowerOfTwo((columns-1)) ) {
 967     AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
 968     return; }
 969   if ( phislices <= 3   )  {
 970     AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
 971     return; }
 972   if  ( phislices > 1000 ) {
 973     AliError("Poisson3D  phislices > 1000 is not allowed (nor wise) ");
 974     return; }
 975
 976   // Solve Poisson's equation in cylindrical coordinates by relaxation technique
 977   // Allow for different size grid spacing in R and Z directions
 978   // Use a binary expansion of the matrix to speed up the solution of the problem
 979
 980   Int_t loops, mplus, mminus, signplus, signminus  ;
 981   Int_t ione = (rows-1)/4 ;
 982   Int_t jone = (columns-1)/4 ;
 983   loops = TMath::Max(ione, jone) ;      // Calculate the number of loops for the binary expansion
 984   loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ;  // Solve for N in 2**N
 985
 986   TMatrixD* arrayofSumChargeDensities[1000] ;    // Create temporary arrays to store low resolution charge arrays
 987
 988   for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
 989
 990   for ( Int_t count = 0 ; count < loops ; count++ ) {      // START the master loop and do the binary expansion
 991
 992     Float_t  tempgridSizeR   =  gridSizeR  * ione ;
 993     Float_t  tempratioPhi    =  ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
 994     Float_t  tempratioZ      =  ratioZ   * ione * ione / ( jone * jone ) ;
 995
 996     std::vector<float> coef1(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
 997     std::vector<float> coef2(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
 998     std::vector<float> coef3(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
 999     std::vector<float> coef4(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1000
1001     for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
1002       Float_t radius = fgkIFCRadius + i*gridSizeR ;
1003       coef1[i] = 1.0 + tempgridSizeR/(2*radius);
1004       coef2[i] = 1.0 - tempgridSizeR/(2*radius);
1005       coef3[i] = tempratioPhi/(radius*radius);
1006       coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
1007     }
1008
1009     for ( Int_t m = 0 ; m < phislices ; m++ ) {
1010       TMatrixD &chargeDensity    = *arrayofChargeDensities[m] ;
1011       TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
1012       for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
1013         Float_t radius = fgkIFCRadius + i*gridSizeR ;
1014         for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
1015           if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1016           else {           // Add up all enclosed charge density contributions within 1/2 unit in all directions
1017             Float_t weight = 0.0 ;
1018             Float_t sum    = 0.0 ;
1019             sumChargeDensity(i,j) = 0.0 ;
1020             for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
1021               for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
1022                 if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
1023                 else
1024                   weight = 1.0 ;
1025                 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1026                 sum += weight*radius ;
1027               }
1028             }
1029             sumChargeDensity(i,j) /= sum ;
1030           }
1031           sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
1032         }
1033       }
1034     }
1035
1036     for ( Int_t k = 1 ; k <= iterations; k++ ) {
1037
1038       // over-relaxation index, >= 1 but < 2
1039       Float_t overRelax   = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1040       Float_t overRelaxM1 = overRelax - 1.0 ;
1041
1042       std::vector<float> overRelaxcoef4(rows) ;  // Do this the standard C++ way to avoid gcc extensions
1043       std::vector<float> overRelaxcoef5(rows) ;  // Do this the standard C++ way to avoid gcc extensions
1044
1045       for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1046         overRelaxcoef4[i] = overRelax * coef4[i] ;
1047         overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
1048       }
1049
1050       for ( Int_t m = 0 ; m < phislices ; m++ ) {
1051
1052         mplus  = m + 1;   signplus  = 1 ;
1053         mminus = m - 1 ;  signminus = 1 ;
1054         if (symmetry==1) {  // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1055           if ( mplus  > phislices-1 ) mplus  = phislices - 2 ;
1056           if ( mminus < 0 )           mminus = 1 ;
1057         }
1058         else if (symmetry==-1) {   // Anti-symmetry in phi
1059           if ( mplus  > phislices-1 ) { mplus  = phislices - 2 ; signplus  = -1 ; }
1060           if ( mminus < 0 )           { mminus = 1 ;             signminus = -1 ; }
1061         }
1062                 else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1063           if ( mplus  > phislices-1 ) mplus  = m + 1 - phislices ;
1064           if ( mminus < 0 )           mminus = m - 1 + phislices ;
1065         }
1066         TMatrixD& arrayV    =  *arrayofArrayV[m] ;
1067         TMatrixD& arrayVP   =  *arrayofArrayV[mplus] ;
1068         TMatrixD& arrayVM   =  *arrayofArrayV[mminus] ;
1069         TMatrixD& sumChargeDensity =  *arrayofSumChargeDensities[m] ;
1070
1071         for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
1072           for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1073
1074             arrayV(i,j) = (   coef2[i]          *   arrayV(i-ione,j)
1075                             + tempratioZ        * ( arrayV(i,j-jone)  +  arrayV(i,j+jone) )
1076                             - overRelaxcoef5[i] *   arrayV(i,j)
1077                             + coef1[i]          *   arrayV(i+ione,j)
1078                             + coef3[i]          * ( signplus*arrayVP(i,j)       +  signminus*arrayVM(i,j) )
1079                             + sumChargeDensity(i,j)
1080                           ) * overRelaxcoef4[i] ;
1081             // Note: over-relax the solution at each step.  This speeds up the convergance.
1082
1083           }
1084         }
1085
1086         if ( k == iterations ) {   // After full solution is achieved, copy low resolution solution into higher res array
1087           for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
1088             for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1089
1090               if ( ione > 1 ) {
1091                 arrayV(i+ione/2,j)                    =  ( arrayV(i+ione,j) + arrayV(i,j)     ) / 2 ;
1092                 if ( i == ione )  arrayV(i-ione/2,j) =  ( arrayV(0,j)       + arrayV(ione,j) ) / 2 ;
1093               }
1094               if ( jone > 1 ) {
1095                 arrayV(i,j+jone/2)                    =  ( arrayV(i,j+jone) + arrayV(i,j) )     / 2 ;
1096                 if ( j == jone )  arrayV(i,j-jone/2) =  ( arrayV(i,0)       + arrayV(i,jone) ) / 2 ;
1097               }
1098               if ( ione > 1 && jone > 1 ) {
1099                 arrayV(i+ione/2,j+jone/2) =  ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
1100                 if ( i == ione ) arrayV(i-ione/2,j-jone/2) =   ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
1101                 if ( j == jone ) arrayV(i-ione/2,j-jone/2) =   ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
1102                 // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
1103               }
1104             }
1105           }
1106         }
1107
1108       }
1109     }
1110
1111     ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
1112     jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
1113
1114   }
1115
1116   //Differentiate V(r) and solve for E(r) using special equations for the first and last row
1117   //Integrate E(r)/E(z) from point of origin to pad plane
1118
1119   for ( Int_t m = 0 ; m < phislices ; m++ ) {
1120     TMatrixD& arrayV    =  *arrayofArrayV[m] ;
1121     TMatrixD& eroverEz  =  *arrayofEroverEz[m] ;
1122
1123     for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {  // Count backwards to facilitate integration over Z
1124
1125       // Differentiate in R
1126       for ( Int_t i = 1 ; i < rows-1 ; i++ )  arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1127       arrayE(0,j)      =  -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1128       arrayE(rows-1,j) =  -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1129       // Integrate over Z
1130       for ( Int_t i = 0 ; i < rows ; i++ ) {
1131         Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
1132         eroverEz(i,j) = 0.0 ;
1133         for ( Int_t k = j ; k < columns ; k++ ) {
1134
1135           eroverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1136           if ( index != 4 )  index = 4; else index = 2 ;
1137         }
1138         if ( index == 4 ) eroverEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1139         if ( index == 2 ) eroverEz(i,j)  +=
1140           (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1141         if ( j == columns-2 ) eroverEz(i,j) =
1142           (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1143         if ( j == columns-1 ) eroverEz(i,j) =  0.0 ;
1144       }
1145     }
1146     // if ( m == 0 ) { TCanvas*  c1 =  new TCanvas("erOverEz","erOverEz",50,50,840,600) ;  c1 -> cd() ;
1147     // eroverEz.Draw("surf") ; } // JT test
1148   }
1149
1150   //Differentiate V(r) and solve for E(phi)
1151   //Integrate E(phi)/E(z) from point of origin to pad plane
1152
1153   for ( Int_t m = 0 ; m < phislices ; m++ ) {
1154
1155     mplus  = m + 1;   signplus  = 1 ;
1156     mminus = m - 1 ;  signminus = 1 ;
1157     if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1158       if ( mplus  > phislices-1 ) mplus  = phislices - 2 ;
1159       if ( mminus < 0 )           mminus = 1 ;
1160     }
1161     else if (symmetry==-1) {       // Anti-symmetry in phi
1162       if ( mplus  > phislices-1 ) { mplus  = phislices - 2 ;  signplus  = -1 ; }
1163       if ( mminus < 0 )           { mminus = 1 ;                    signminus = -1 ; }
1164     }
1165     else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
1166       if ( mplus  > phislices-1 ) mplus  = m + 1 - phislices ;
1167       if ( mminus < 0 )           mminus = m - 1 + phislices ;
1168     }
1169     TMatrixD &arrayVP     =  *arrayofArrayV[mplus] ;
1170     TMatrixD &arrayVM     =  *arrayofArrayV[mminus] ;
1171     TMatrixD &ePhioverEz  =  *arrayofEPhioverEz[m] ;
1172     for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1173       // Differentiate in Phi
1174       for ( Int_t i = 0 ; i < rows ; i++ ) {
1175         Float_t radius = fgkIFCRadius + i*gridSizeR ;
1176         arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
1177       }
1178       // Integrate over Z
1179       for ( Int_t i = 0 ; i < rows ; i++ ) {
1180         Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
1181         ePhioverEz(i,j) = 0.0 ;
1182         for ( Int_t k = j ; k < columns ; k++ ) {
1183
1184           ePhioverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1185           if ( index != 4 )  index = 4; else index = 2 ;
1186         }
1187         if ( index == 4 ) ePhioverEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1188         if ( index == 2 ) ePhioverEz(i,j)  +=
1189           (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1190         if ( j == columns-2 ) ePhioverEz(i,j) =
1191           (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1192         if ( j == columns-1 ) ePhioverEz(i,j) =  0.0 ;
1193       }
1194     }
1195     // if ( m == 5 ) { TCanvas* c2 =  new TCanvas("arrayE","arrayE",50,50,840,600) ;  c2 -> cd() ;
1196     // arrayE.Draw("surf") ; } // JT test
1197   }
1198
1199
1200   // Differentiate V(r) and solve for E(z) using special equations for the first and last row
1201   // Integrate (E(z)-Ezstd) from point of origin to pad plane
1202
1203   for ( Int_t m = 0 ; m < phislices ; m++ ) {
1204     TMatrixD& arrayV   =  *arrayofArrayV[m] ;
1205     TMatrixD& deltaEz  =  *arrayofDeltaEz[m] ;
1206
1207     // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1208     for ( Int_t i = 0 ; i < rows ; i++) {
1209       for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1210       arrayE(i,0)         =  -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1211       arrayE(i,columns-1) =  -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1212     }
1213
1214     for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {  // Count backwards to facilitate integration over Z
1215       // Integrate over Z
1216       for ( Int_t i = 0 ; i < rows ; i++ ) {
1217         Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
1218         deltaEz(i,j) = 0.0 ;
1219         for ( Int_t k = j ; k < columns ; k++ ) {
1220           deltaEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k) ;
1221           if ( index != 4 )  index = 4; else index = 2 ;
1222         }
1223         if ( index == 4 ) deltaEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1) ;
1224         if ( index == 2 ) deltaEz(i,j)  +=
1225           (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
1226         if ( j == columns-2 ) deltaEz(i,j) =
1227           (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
1228         if ( j == columns-1 ) deltaEz(i,j) =  0.0 ;
1229       }
1230     }
1231     // if ( m == 0 ) { TCanvas*  c1 =  new TCanvas("erOverEz","erOverEz",50,50,840,600) ;  c1 -> cd() ;
1232     // eroverEz.Draw("surf") ; } // JT test
1233
1234     // calculate z distortion from the integrated Delta Ez residuals
1235     // and include the aquivalence (Volt to cm) of the ROC shift !!
1236
1237     for ( Int_t j = 0 ; j < columns ; j++ )  {
1238       for ( Int_t i = 0 ; i < rows ; i++ ) {
1239
1240         // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1241         deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
1242
1243         // ROC Potential in cm aquivalent
1244         Double_t dzROCShift =  arrayV(i, columns -1)/ezField;
1245         if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift;  // add the ROC misaligment
1246
1247       }
1248     }
1249
1250   } // end loop over phi
1251
1252
1253
1254   for ( Int_t k = 0 ; k < phislices ; k++ )
1255     {
1256       arrayofSumChargeDensities[k]->Delete() ;
1257     }
1258
1259
1260
1261   arrayE.Clear();
1262 }
1263
1264
1265 Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1266   //
1267   // Helperfunction: Check if integer is a power of 2
1268   //
1269   Int_t j = 0;
1270   while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
1271   if ( j == 1 ) return(1) ;  // True
1272   return(0) ;                // False
1273 }
1274
1275
1276 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
1277   //
1278   // Fit the track parameters - without and with distortion
1279   // 1. Space points in the TPC are simulated along the trajectory
1280   // 2. Space points distorted
1281   // 3. Fits  the non distorted and distroted track to the reference plane at refX
1282   // 4. For visualization and debugging  purposes the space points and tracks can be stored  in the tree - using the TTreeSRedirector functionality
1283   //
1284   // trackIn   - input track parameters
1285   // refX     - reference X to fit the track
1286   // dir      - direction - out=1 or in=-1
1287   // pcstream -  debug streamer to check the results
1288   //
1289   // see AliExternalTrackParam.h documentation:
1290   // track1.fP[0] - local y (rphi)
1291   // track1.fP[1] - z
1292   // track1.fP[2] - sinus of local inclination angle
1293   // track1.fP[3] - tangent of deep angle
1294   // track1.fP[4] - 1/pt
1295
1296   AliTPCROC * roc = AliTPCROC::Instance();
1297   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
1298   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
1299   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1300   const Double_t kMaxSnp = 0.85;
1301   const Double_t kSigmaY=0.1;
1302   const Double_t kSigmaZ=0.1;
1303   const Double_t kMaxR=500;
1304   const Double_t kMaxZ=500;
1305   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1306   Int_t npoints1=0;
1307   Int_t npoints2=0;
1308
1309   AliExternalTrackParam  track(trackIn); //
1310   // generate points
1311   AliTrackPointArray pointArray0(npoints0);
1312   AliTrackPointArray pointArray1(npoints0);
1313   Double_t xyz[3];
1314   if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
1315   //
1316   // simulate the track
1317   Int_t npoints=0;
1318   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
1319   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
1320     if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
1321     track.GetXYZ(xyz);
1322     xyz[0]+=gRandom->Gaus(0,0.00005);
1323     xyz[1]+=gRandom->Gaus(0,0.00005);
1324     xyz[2]+=gRandom->Gaus(0,0.00005);
1325     if (TMath::Abs(track.GetZ())>kMaxZ) break;
1326     if (TMath::Abs(track.GetX())>kMaxR) break;
1327     AliTrackPoint pIn0;                               // space point
1328     AliTrackPoint pIn1;
1329     Int_t sector= (xyz[2]>0)? 0:18;
1330     pointArray0.GetPoint(pIn0,npoints);
1331     pointArray1.GetPoint(pIn1,npoints);
1332     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
1333     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
1334     DistortPoint(distPoint, sector);
1335     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
1336     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
1337     //
1338     track.Rotate(alpha);
1339     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
1340     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
1341     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
1342     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
1343     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
1344     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
1345     pointArray0.AddPoint(npoints, &pIn0);
1346     pointArray1.AddPoint(npoints, &pIn1);
1347     npoints++;
1348     if (npoints>=npoints0) break;
1349   }
1350   if (npoints<npoints0/2) return 0;
1351   //
1352   // refit track
1353   //
1354   AliExternalTrackParam *track0=0;
1355   AliExternalTrackParam *track1=0;
1356   AliTrackPoint   point1,point2,point3;
1357   if (dir==1) {  //make seed inner
1358     pointArray0.GetPoint(point1,1);
1359     pointArray0.GetPoint(point2,30);
1360     pointArray0.GetPoint(point3,60);
1361   }
1362   if (dir==-1){ //make seed outer
1363     pointArray0.GetPoint(point1,npoints-60);
1364     pointArray0.GetPoint(point2,npoints-30);
1365     pointArray0.GetPoint(point3,npoints-1);
1366   }
1367   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
1368   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
1369
1370   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
1371     Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
1372     //
1373     AliTrackPoint pIn0;
1374     AliTrackPoint pIn1;
1375     pointArray0.GetPoint(pIn0,ipoint);
1376     pointArray1.GetPoint(pIn1,ipoint);
1377     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
1378     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
1379     //
1380     if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1381     if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1382     if (TMath::Abs(track0->GetZ())>kMaxZ) break;
1383     if (TMath::Abs(track0->GetX())>kMaxR) break;
1384     if (TMath::Abs(track1->GetZ())>kMaxZ) break;
1385     if (TMath::Abs(track1->GetX())>kMaxR) break;
1386
1387     track.GetXYZ(xyz);  // distorted track also propagated to the same reference radius
1388     //
1389     Double_t pointPos[2]={0,0};
1390     Double_t pointCov[3]={0,0,0};
1391     pointPos[0]=prot0.GetY();//local y
1392     pointPos[1]=prot0.GetZ();//local z
1393     pointCov[0]=prot0.GetCov()[3];//simay^2
1394     pointCov[1]=prot0.GetCov()[4];//sigmayz
1395     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
1396     if (!track0->Update(pointPos,pointCov)) break;
1397     //
1398     Double_t deltaX=prot1.GetX()-prot0.GetX();   // delta X
1399     Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp()));  // deltaY due  delta X
1400     Double_t deltaZX=deltaX*track1->GetTgl();                           // deltaZ due  delta X
1401
1402     pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
1403     pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
1404     pointCov[0]=prot1.GetCov()[3];//simay^2
1405     pointCov[1]=prot1.GetCov()[4];//sigmayz
1406     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
1407     if (!track1->Update(pointPos,pointCov)) break;
1408     npoints1++;
1409     npoints2++;
1410   }
1411   if (npoints2<npoints)  return 0;
1412   AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
1413   track1->Rotate(track0->GetAlpha());
1414   AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
1415
1416   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
1417     "point0.="<<&pointArray0<<   //  points
1418     "point1.="<<&pointArray1<<   //  distorted points
1419     "trackIn.="<<&track<<       //  original track
1420     "track0.="<<track0<<         //  fitted track
1421     "track1.="<<track1<<         //  fitted distorted track
1422     "\n";
1423   new(&trackIn) AliExternalTrackParam(*track0);
1424   delete track0;
1425   return track1;
1426 }
1427
1428
1429
1430
1431
1432 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
1433   //
1434   // create the distortion tree on a mesh with granularity given by step
1435   // return the tree with distortions at given position
1436   // Map is created on the mesh with given step size
1437   //
1438   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
1439   Float_t xyz[3];
1440   for (Double_t x= -250; x<250; x+=step){
1441     for (Double_t y= -250; y<250; y+=step){
1442       Double_t r    = TMath::Sqrt(x*x+y*y);
1443       if (r<80) continue;
1444       if (r>250) continue;
1445       for (Double_t z= -250; z<250; z+=step){
1446         Int_t roc=(z>0)?0:18;
1447         xyz[0]=x;
1448         xyz[1]=y;
1449         xyz[2]=z;
1450         Double_t phi  = TMath::ATan2(y,x);
1451         DistortPoint(xyz,roc);
1452         Double_t r1    = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1453         Double_t phi1  = TMath::ATan2(xyz[1],xyz[0]);
1454         if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
1455         if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
1456         Double_t dx = xyz[0]-x;
1457         Double_t dy = xyz[1]-y;
1458         Double_t dz = xyz[2]-z;
1459         Double_t dr=r1-r;
1460         Double_t drphi=(phi1-phi)*r;
1461         (*pcstream)<<"distortion"<<
1462           "x="<<x<<           // original position
1463           "y="<<y<<
1464           "z="<<z<<
1465           "r="<<r<<
1466           "phi="<<phi<<
1467           "x1="<<xyz[0]<<      // distorted position
1468           "y1="<<xyz[1]<<
1469           "z1="<<xyz[2]<<
1470           "r1="<<r1<<
1471           "phi1="<<phi1<<
1472           //
1473           "dx="<<dx<<          // delta position
1474           "dy="<<dy<<
1475           "dz="<<dz<<
1476           "dr="<<dr<<
1477           "drphi="<<drphi<<
1478           "\n";
1479       }
1480     }
1481   }
1482   delete pcstream;
1483   TFile f(Form("correction%s.root",GetName()));
1484   TTree * tree = (TTree*)f.Get("distortion");
1485   TTree * tree2= tree->CopyTree("1");
1486   tree2->SetName(Form("dist%s",GetName()));
1487   tree2->SetDirectory(0);
1488   delete tree;
1489   return tree2;
1490 }
1491
1492
1493
1494
1495 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
1496   //
1497   // Make a fit tree:
1498   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1499   // calculates partial distortions
1500   // Partial distortion is stored in the resulting tree
1501   // Output is storred in the file distortion_<dettype>_<partype>.root
1502   // Partial  distortion is stored with the name given by correction name
1503   //
1504   //
1505   // Parameters of function:
1506   // input     - input tree
1507   // dtype     - distortion type 0 - ITSTPC,  1 -TPCTRD, 2 - TPCvertex
1508   // ppype     - parameter type
1509   // corrArray - array with partial corrections
1510   // step      - skipe entries  - if 1 all entries processed - it is slow
1511   // debug     0 if debug on also space points dumped - it is slow
1512
1513   const Double_t kMaxSnp = 0.85;
1514   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1515   //  const Double_t kB2C=-0.299792458e-3;
1516   const Int_t kMinEntries=50;
1517   Double_t phi,theta, snp, mean,rms, entries;
1518   tinput->SetBranchAddress("theta",&theta);
1519   tinput->SetBranchAddress("phi", &phi);
1520   tinput->SetBranchAddress("snp",&snp);
1521   tinput->SetBranchAddress("mean",&mean);
1522   tinput->SetBranchAddress("rms",&rms);
1523   tinput->SetBranchAddress("entries",&entries);
1524   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
1525   //
1526   Int_t nentries=tinput->GetEntries();
1527   Int_t ncorr=corrArray->GetEntries();
1528   Double_t corrections[100]={0}; //
1529   Double_t tPar[5];
1530   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1531   Double_t refX=0;
1532   Int_t dir=0;
1533   if (dtype==0) {refX=85.; dir=-1;}
1534   if (dtype==1) {refX=275.; dir=1;}
1535   if (dtype==2) {refX=85.; dir=-1;}
1536   if (dtype==3) {refX=360.; dir=-1;}
1537   //
1538   for (Int_t ientry=0; ientry<nentries; ientry+=step){
1539     tinput->GetEntry(ientry);
1540     if (TMath::Abs(snp)>kMaxSnp) continue;
1541     tPar[0]=0;
1542     tPar[1]=theta*refX;
1543     tPar[2]=snp;
1544     tPar[3]=theta;
1545     tPar[4]=(gRandom->Rndm()-0.5)*0.02;  // should be calculated - non equal to 0
1546     Double_t bz=AliTrackerBase::GetBz();
1547     if (refX>10. && TMath::Abs(bz)>0.1 )  tPar[4]=snp/(refX*bz*kB2C*2);
1548     tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
1549     AliExternalTrackParam track(refX,phi,tPar,cov);
1550     Double_t xyz[3];
1551     track.GetXYZ(xyz);
1552     Int_t id=0;
1553     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1554     if (ptype==4 &&bz<0) mean*=-1;  // interpret as curvature
1555     (*pcstream)<<"fit"<<
1556       "bz="<<bz<<         // magnetic filed used
1557       "dtype="<<dtype<<   // detector match type
1558       "ptype="<<ptype<<   // parameter type
1559       "theta="<<theta<<   // theta
1560       "phi="<<phi<<       // phi
1561       "snp="<<snp<<       // snp
1562       "mean="<<mean<<     // mean dist value
1563       "rms="<<rms<<       // rms
1564       "gx="<<xyz[0]<<         // global position at reference
1565       "gy="<<xyz[1]<<         // global position at reference
1566       "gz="<<xyz[2]<<         // global position at reference
1567       "dRrec="<<dRrec<<      // delta Radius in reconstruction
1568       "id="<<id<<             // track id
1569       "entries="<<entries;// number of entries in bin
1570     //
1571     for (Int_t icorr=0; icorr<ncorr; icorr++) {
1572       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1573       corrections[icorr]=0;
1574       if (entries>kMinEntries){
1575         AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1576         AliExternalTrackParam *trackOut = 0;
1577         if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1578         if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1579         if (dtype==0) {refX=85.; dir=-1;}
1580         if (dtype==1) {refX=275.; dir=1;}
1581         if (dtype==2) {refX=0; dir=-1;}
1582         if (dtype==3) {refX=360.; dir=-1;}
1583         //
1584         if (trackOut){
1585           AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1586           trackOut->Rotate(trackIn.GetAlpha());
1587           trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1588           //
1589           corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1590           delete trackOut;
1591         }else{
1592           corrections[icorr]=0;
1593         }
1594         if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature
1595       }
1596       Double_t dRdummy=0;
1597       (*pcstream)<<"fit"<<
1598         Form("%s=",corr->GetName())<<corrections[icorr]<<   // dump correction value
1599         Form("dR%s=",corr->GetName())<<dRdummy;   // dump dummy correction value not needed for tracks
1600                                                   // for points it is neccessary
1601     }
1602     (*pcstream)<<"fit"<<"\n";
1603   }
1604   delete pcstream;
1605 }
1606
1607
1608
1609 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1610   //
1611   // Make a laser fit tree for global minimization
1612   //
1613   const Double_t cutErrY=0.1;
1614   const Double_t cutErrZ=0.1;
1615   const Double_t kEpsilon=0.00000001;
1616   TVectorD *vecdY=0;
1617   TVectorD *vecdZ=0;
1618   TVectorD *veceY=0;
1619   TVectorD *veceZ=0;
1620   AliTPCLaserTrack *ltr=0;
1621   AliTPCLaserTrack::LoadTracks();
1622   tree->SetBranchAddress("dY.",&vecdY);
1623   tree->SetBranchAddress("dZ.",&vecdZ);
1624   tree->SetBranchAddress("eY.",&veceY);
1625   tree->SetBranchAddress("eZ.",&veceZ);
1626   tree->SetBranchAddress("LTr.",&ltr);
1627   Int_t entries= tree->GetEntries();
1628   TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1629   Double_t bz=AliTrackerBase::GetBz();
1630   //
1631
1632   for (Int_t ientry=0; ientry<entries; ientry++){
1633     tree->GetEntry(ientry);
1634     if (!ltr->GetVecGX()){
1635       ltr->UpdatePoints();
1636     }
1637     TVectorD * delta= (itype==0)? vecdY:vecdZ;
1638     TVectorD * err= (itype==0)? veceY:veceZ;
1639
1640     for (Int_t irow=0; irow<159; irow++){
1641       Int_t nentries = 1000;
1642       if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1643       if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1644       Int_t dtype=4;
1645       Double_t phi   =(*ltr->GetVecPhi())[irow];
1646       Double_t theta =ltr->GetTgl();
1647       Double_t mean=delta->GetMatrixArray()[irow];
1648       Double_t gx=0,gy=0,gz=0;
1649       Double_t snp = (*ltr->GetVecP2())[irow];
1650       Double_t rms = 0.1+err->GetMatrixArray()[irow];
1651       gx = (*ltr->GetVecGX())[irow];
1652       gy = (*ltr->GetVecGY())[irow];
1653       gz = (*ltr->GetVecGZ())[irow];
1654       Int_t bundle= ltr->GetBundle();
1655       Double_t dRrec=0;
1656       //
1657       // get delta R used in reconstruction
1658       AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
1659       AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1660       const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1661       Double_t xyz0[3]={gx,gy,gz};
1662       Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1663       //
1664       // old ExB correction
1665       //
1666       if(recoParam&&recoParam->GetUseExBCorrection()) {
1667         Double_t xyz1[3]={gx,gy,gz};
1668         calib->GetExB()->Correct(xyz0,xyz1);
1669         Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1670         dRrec=oldR-newR;
1671       }
1672       if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1673         Float_t xyz1[3]={gx,gy,gz};
1674         Int_t sector=(gz>0)?0:18;
1675         correction->CorrectPoint(xyz1, sector);
1676         Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1677         dRrec=oldR-newR;
1678       }
1679
1680
1681       (*pcstream)<<"fit"<<
1682         "bz="<<bz<<         // magnetic filed used
1683         "dtype="<<dtype<<   // detector match type
1684         "ptype="<<itype<<   // parameter type
1685         "theta="<<theta<<   // theta
1686         "phi="<<phi<<       // phi
1687         "snp="<<snp<<       // snp
1688         "mean="<<mean<<     // mean dist value
1689         "rms="<<rms<<       // rms
1690         "gx="<<gx<<         // global position
1691         "gy="<<gy<<         // global position
1692         "gz="<<gz<<         // global position
1693         "dRrec="<<dRrec<<      // delta Radius in reconstruction
1694         "id="<<bundle<<     //bundle
1695         "entries="<<nentries;// number of entries in bin
1696       //
1697       //
1698       Double_t ky = TMath::Tan(TMath::ASin(snp));
1699       Int_t ncorr = corrArray->GetEntries();
1700       Double_t r0   = TMath::Sqrt(gx*gx+gy*gy);
1701       Double_t phi0 = TMath::ATan2(gy,gx);
1702       Double_t distortions[1000]={0};
1703       Double_t distortionsR[1000]={0};
1704       for (Int_t icorr=0; icorr<ncorr; icorr++) {
1705         AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1706         Float_t distPoint[3]={gx,gy,gz};
1707         Int_t sector= (gz>0)? 0:18;
1708         if (r0>80){
1709           corr->DistortPoint(distPoint, sector);
1710         }
1711         // Double_t value=distPoint[2]-gz;
1712         if (itype==0){
1713           Double_t r1   = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1714           Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1715           Double_t drphi= r0*(phi1-phi0);
1716           Double_t dr   = r1-r0;
1717           distortions[icorr]  = drphi-ky*dr;
1718           distortionsR[icorr] = dr;
1719         }
1720         (*pcstream)<<"fit"<<
1721           Form("%s=",corr->GetName())<<distortions[icorr]<<    // dump correction value
1722           Form("dR%s=",corr->GetName())<<distortionsR[icorr];   // dump correction R  value
1723       }
1724       (*pcstream)<<"fit"<<"\n";
1725     }
1726   }
1727   delete pcstream;
1728 }
1729
1730
1731
1732 void   AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
1733   //
1734   // make a distortion map out ou fthe residual histogram
1735   // Results are written to the debug streamer - pcstream
1736   // Parameters:
1737   //   his0       - input (4D) residual histogram
1738   //   pcstream   - file to write the tree
1739   //   run        - run number
1740   // marian.ivanov@cern.ch
1741   const Int_t kMinEntries=50;
1742   Int_t nbins1=his0->GetAxis(1)->GetNbins();
1743   Int_t first1=his0->GetAxis(1)->GetFirst();
1744   Int_t last1 =his0->GetAxis(1)->GetLast();
1745   //
1746   Double_t bz=AliTrackerBase::GetBz();
1747   Int_t idim[4]={0,1,2,3};
1748   for (Int_t ibin1=first1; ibin1<last1; ibin1++){   //axis 1 - theta
1749     //
1750     his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1751     Double_t       x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1752     THnSparse * his1 = his0->Projection(4,idim);  // projected histogram according range1
1753     Int_t nbins3     = his1->GetAxis(3)->GetNbins();
1754     Int_t first3     = his1->GetAxis(3)->GetFirst();
1755     Int_t last3      = his1->GetAxis(3)->GetLast();
1756     //
1757
1758     for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){   // axis 3 - local angle
1759       his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1760       Double_t      x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1761       if (ibin3<first3) {
1762         his1->GetAxis(3)->SetRangeUser(-1,1);
1763         x3=0;
1764       }
1765       THnSparse * his3= his1->Projection(4,idim);         //projected histogram according selection 3
1766       Int_t nbins2    = his3->GetAxis(2)->GetNbins();
1767       Int_t first2    = his3->GetAxis(2)->GetFirst();
1768       Int_t last2     = his3->GetAxis(2)->GetLast();
1769       //
1770       for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1771         his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1772         Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1773         TH1 * hisDelta = his3->Projection(0);
1774         //
1775         Double_t entries = hisDelta->GetEntries();
1776         Double_t mean=0, rms=0;
1777         if (entries>kMinEntries){
1778           mean    = hisDelta->GetMean();
1779           rms = hisDelta->GetRMS();
1780         }
1781         (*pcstream)<<hname<<
1782           "run="<<run<<
1783           "bz="<<bz<<
1784           "theta="<<x1<<
1785           "phi="<<x2<<
1786           "snp="<<x3<<
1787           "entries="<<entries<<
1788           "mean="<<mean<<
1789           "rms="<<rms<<
1790           "\n";
1791         delete hisDelta;
1792         printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1793       }
1794       delete his3;
1795     }
1796     delete his1;
1797   }
1798 }
1799
1800
1801
1802
1803
1804 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1805   //
1806   // Store object in the OCDB
1807   // By default the object is stored in the current directory
1808   // default comment consit of user name and the date
1809   //
1810   TString ocdbStorage="";
1811   ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1812   AliCDBMetaData *metaData= new AliCDBMetaData();
1813   metaData->SetObjectClassName("AliTPCCorrection");
1814   metaData->SetResponsible("Marian Ivanov");
1815   metaData->SetBeamPeriod(1);
1816   metaData->SetAliRootVersion("05-25-01"); //root version
1817   TString userName=gSystem->GetFromPipe("echo $USER");
1818   TString date=gSystem->GetFromPipe("date");
1819
1820   if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1821   if (comment) metaData->SetComment(comment);
1822   AliCDBId* id1=NULL;
1823   id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1824   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1825   gStorage->Put(this, (*id1), metaData);
1826 }
1827
1828
1829 void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
1830   //
1831   // Fast method to simulate the influence of the given distortion on the vertex reconstruction
1832   //
1833
1834   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
1835   if (!magF) AliError("Magneticd field - not initialized");
1836   Double_t bz = magF->SolenoidField(); //field in kGauss
1837   printf("bz: %lf\n",bz);
1838   AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
1839
1840   TObjArray   aTrk;              // Original Track array of Aside
1841   TObjArray   daTrk;             // Distorted Track array of A side
1842   UShort_t    *aId = new UShort_t[nTracks];      // A side Track ID
1843   TObjArray   cTrk;
1844   TObjArray   dcTrk;
1845   UShort_t    *cId = new UShort_t [nTracks];
1846   Int_t id=0;
1847   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1848   TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
1849   fpt.SetParameters(7.24,0.120);
1850   fpt.SetNpx(10000);
1851   for(Int_t nt=0; nt<nTracks; nt++){
1852     Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
1853     Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
1854     Double_t pt = fpt.GetRandom(); // momentum for f1
1855     //   printf("phi %lf  eta %lf pt %lf\n",phi,eta,pt);
1856     Short_t sign=1;
1857     if(gRandom->Rndm() < 0.5){
1858       sign =1;
1859     }else{
1860       sign=-1;
1861     }
1862
1863     Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1864     Double_t pxyz[3];
1865     pxyz[0]=pt*TMath::Cos(phi);
1866     pxyz[1]=pt*TMath::Sin(phi);
1867     pxyz[2]=pt*TMath::Tan(theta);
1868     Double_t cv[21]={0};
1869     AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
1870
1871     Double_t refX=1.;
1872     Int_t dir=-1;
1873     AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir,  NULL);
1874     if (!td) continue;
1875     if (pcstream) (*pcstream)<<"track"<<
1876       "eta="<<eta<<
1877       "theta="<<theta<<
1878       "tOrig.="<<t<<
1879       "td.="<<td<<
1880       "\n";
1881     if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
1882       if (td){
1883         daTrk.AddLast(td);
1884         aTrk.AddLast(t);
1885         Int_t nn=aTrk.GetEntriesFast();
1886         aId[nn]=id;
1887       }
1888     }else if(( eta<-0.07 )&&( eta>-etaCuts )){
1889       if (td){
1890         dcTrk.AddLast(td);
1891         cTrk.AddLast(t);
1892         Int_t nn=cTrk.GetEntriesFast();
1893         cId[nn]=id;
1894       }
1895     }
1896     id++;
1897   }// end of track loop
1898
1899   vertexer->SetTPCMode();
1900   vertexer->SetConstraintOff();
1901
1902   aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));
1903   avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
1904   cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));
1905   cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
1906   if (pcstream) (*pcstream)<<"vertex"<<
1907     "x="<<orgVertex[0]<<
1908     "y="<<orgVertex[1]<<
1909     "z="<<orgVertex[2]<<
1910     "av.="<<&aV<<              // distorted vertex A side
1911     "cv.="<<&cV<<              // distroted vertex C side
1912     "avO.="<<&avOrg<<         // original vertex A side
1913     "cvO.="<<&cvOrg<<
1914     "\n";
1915   delete []aId;
1916   delete []cId;
1917 }
1918
1919 void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
1920   //
1921   // make correction available for visualization using
1922   // TFormula, TFX and TTree::Draw
1923   // important in order to check corrections and also compute dervied variables
1924   // e.g correction partial derivatives
1925   //
1926   // NOTE - class is not owner of correction
1927   //
1928   if (!fgVisualCorrection) fgVisualCorrection=new TObjArray;
1929   if (position!=0&&position>=fgVisualCorrection->GetEntriesFast())
1930     fgVisualCorrection->Expand(position*2);
1931   fgVisualCorrection->AddAt(corr, position);
1932 }
1933
1934
1935
1936 Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t r, Double_t kZ, Int_t axisType, Int_t corrType){
1937   //
1938   // calculate the correction at given position - check the geffCorr
1939   //
1940   if (!fgVisualCorrection) return 0;
1941   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1942   if (!corr) return 0;
1943
1944   Double_t phi=sector*TMath::Pi()/9.;
1945   Double_t gx = r*TMath::Cos(phi);
1946   Double_t gy = r*TMath::Sin(phi);
1947   Double_t gz = r*kZ;
1948   Int_t nsector=(gz>0) ? 0:18;
1949   //
1950   //
1951   //
1952   Float_t distPoint[3]={gx,gy,gz};
1953   corr->DistortPoint(distPoint, nsector);
1954   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1955   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1956   Double_t phi0=TMath::ATan2(gy,gx);
1957   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1958   if (axisType==0) return r1-r0;
1959   if (axisType==1) return (phi1-phi0)*r0;
1960   if (axisType==2) return distPoint[2]-gz;
1961   return phi1-phi0;
1962 }
1963
1964 Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
1965   //
1966   // return correction at given x,y,z
1967   //
1968   if (!fgVisualCorrection) return 0;
1969   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1970   if (!corr) return 0;
1971   Double_t phi0= TMath::ATan2(gy,gx);
1972   Int_t nsector=(gz>0) ? 0:18;
1973   Float_t distPoint[3]={gx,gy,gz};
1974   corr->DistortPoint(distPoint, nsector);
1975   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1976   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1977   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1978   if (axisType==0) return r1-r0;
1979   if (axisType==1) return (phi1-phi0)*r0;
1980   if (axisType==2) return distPoint[2]-gz;
1981   return phi1-phi0;
1982 }