Two minor bugfixes
[u/mrichter/AliRoot.git] / TPC / AliTPCCorrection.cxx
index 3fcc1b5..547d24f 100644 (file)
 #include <AliCDBStorage.h>
 #include <AliCDBId.h>
 #include <AliCDBMetaData.h>
+#include "TVectorD.h"
+#include "AliTPCParamSR.h"
 
-#include  "TRandom.h"
-#include  "AliExternalTrackParam.h"
-#include  "AliTrackPointArray.h"
-#include  "TDatabasePDG.h"
-#include  "AliTrackerBase.h"
-#include  "AliTPCROC.h"
-#include  "THnSparse.h"
+#include "AliTPCCorrection.h"
+#include "AliLog.h"
 
+#include "AliExternalTrackParam.h"
+#include "AliTrackPointArray.h"
+#include "TDatabasePDG.h"
+#include "AliTrackerBase.h"
+#include "AliTPCROC.h"
+#include "THnSparse.h"
+
+#include "AliTPCLaserTrack.h"
+#include "AliESDVertex.h"
+#include "AliVertexerTracks.h"
+#include "TDatabasePDG.h"
+#include "TF1.h"
+#include "TRandom.h"
+
+#include "TDatabasePDG.h"
+
+#include "AliTPCTransform.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCExB.h"
+
+#include "AliTPCRecoParam.h"
 
-#include "AliTPCCorrection.h"
 
 ClassImp(AliTPCCorrection)
 
+
+TObjArray *AliTPCCorrection::fgVisualCorrection=0;
+// instance of correction for visualization
+
+
 // FIXME: the following values should come from the database
-const Double_t AliTPCCorrection::fgkTPC_Z0   =249.7;     // nominal gating grid position 
-const Double_t AliTPCCorrection::fgkIFCRadius= 83.06;    // Mean Radius of the Inner Field Cage ( 82.43 min,  83.70 max) (cm)
-const Double_t AliTPCCorrection::fgkOFCRadius=254.5;     // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
-const Double_t AliTPCCorrection::fgkZOffSet  = 0.2;      // Offset from CE: calculate all distortions closer to CE as if at this point
-const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
-const Double_t AliTPCCorrection::fgkGG       =-70.0;     // Gating Grid voltage (volts)
+const Double_t AliTPCCorrection::fgkTPCZ0    = 249.7;     // nominal gating grid position 
+const Double_t AliTPCCorrection::fgkIFCRadius=  83.5;     // radius which renders the "18 rod manifold" best -> compare calc. of Jim Thomas
+// compare gkIFCRadius=  83.05: Mean Radius of the Inner Field Cage ( 82.43 min,  83.70 max) (cm)
+const Double_t AliTPCCorrection::fgkOFCRadius= 254.5;     // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
+const Double_t AliTPCCorrection::fgkZOffSet  =   0.2;     // Offset from CE: calculate all distortions closer to CE as if at this point
+const Double_t AliTPCCorrection::fgkCathodeV = -100000.0; // Cathode Voltage (volts)
+const Double_t AliTPCCorrection::fgkGG       =     -70.0; // Gating Grid voltage (volts)
+
+const Double_t  AliTPCCorrection::fgkdvdE = 0.0024; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
 
+const Double_t AliTPCCorrection::fgkEM = -1.602176487e-19/9.10938215e-31; // charge/mass in [C/kg]
+const Double_t AliTPCCorrection::fgke0 = 8.854187817e-12;                 // vacuum permittivity [A·s/(V·m)]
 
 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] =  {   
-84.0,   84.5,   85.0,   85.5,   86.0,  87.0,    88.0,
-90.0,   92.0,   94.0,   96.0,   98.0,  100.0,  102.0,  104.0,  106.0,  108.0, 
-110.0,  112.0,  114.0,  116.0,  118.0,  120.0,  122.0,  124.0,  126.0,  128.0, 
-130.0,  132.0,  134.0,  136.0,  138.0,  140.0,  142.0,  144.0,  146.0,  148.0, 
-150.0,  152.0,  154.0,  156.0,  158.0,  160.0,  162.0,  164.0,  166.0,  168.0, 
-170.0,  172.0,  174.0,  176.0,  178.0,  180.0,  182.0,  184.0,  186.0,  188.0,
-190.0,  192.0,  194.0,  196.0,  198.0,  200.0,  202.0,  204.0,  206.0,  208.0,
-210.0,  212.0,  214.0,  216.0,  218.0,  220.0,  222.0,  224.0,  226.0,  228.0,
-230.0,  232.0,  234.0,  236.0,  238.0,  240.0,  242.0,  244.0,  246.0,  248.0,
-249.0,  249.5,  250.0,  251.5,  252.0  } ;
-  
+  83.06,   83.5,   84.0,   84.5,   85.0,   85.5,   86.0,   86.5,      
+   87.0,   87.5,   88.0,   88.5,   89.0,   89.5,   90.0,   90.5,   91.0,   92.0,
+   93.0,   94.0,   95.0,   96.0,   98.0,  100.0,  102.0,  104.0,  106.0,  108.0, 
+  110.0,  112.0,  114.0,  116.0,  118.0,  120.0,  122.0,  124.0,  126.0,  128.0, 
+  130.0,  132.0,  134.0,  136.0,  138.0,  140.0,  142.0,  144.0,  146.0,  148.0, 
+  150.0,  152.0,  154.0,  156.0,  158.0,  160.0,  162.0,  164.0,  166.0,  168.0, 
+  170.0,  172.0,  174.0,  176.0,  178.0,  180.0,  182.0,  184.0,  186.0,  188.0,
+  190.0,  192.0,  194.0,  196.0,  198.0,  200.0,  202.0,  204.0,  206.0,  208.0,
+  210.0,  212.0,  214.0,  216.0,  218.0,  220.0,  222.0,  224.0,  226.0,  228.0,
+  230.0,  232.0,  234.0,  236.0,  238.0,  240.0,  241.0,  242.0,  243.0,  244.0,  
+  245.0,  245.5,  246.0,  246.5,  247.0,  247.5,  248.0,  248.5,  249.0,  249.5,  
+  250.0,  250.5,  251.0,  251.5,  252.0,  252.5,  253.0,  253.5,  254.0,  254.5 
+} ;
 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ]     =   { 
 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
@@ -118,19 +148,33 @@ const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ]     =   {
 
 
 AliTPCCorrection::AliTPCCorrection() 
-  : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
+  : TNamed("correction_unity","unity"),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
 {
   //
   // default constructor
   //
+  if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
+
+  // Initialization of interpolation points 
+  for (Int_t i = 0; i<kNPhi; i++) {
+    fgkPhiList[i] = TMath::TwoPi()*i/(kNPhi-1);    
+  }
+
 }
 
 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
-: TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
+: TNamed(name,title),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
 {
   //
   // default constructor, that set the name and title
   //
+  if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
+
+  // Initialization of interpolation points 
+  for (Int_t i = 0; i<kNPhi; i++) {
+    fgkPhiList[i] = TMath::TwoPi()*i/(kNPhi-1);    
+  }
+
 }
 
 AliTPCCorrection::~AliTPCCorrection() {
@@ -240,7 +284,8 @@ TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
   // in respect to position z within the XY plane.
   // The histogramm has nx times ny entries. 
   //
-  
+  AliTPCParam* tpcparam = new AliTPCParamSR;
+
   TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
                     nx,-250.,250.,ny,-250.,250.);
   Float_t x[3],dx[3];
@@ -252,7 +297,7 @@ TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
       x[0]=h->GetXaxis()->GetBinCenter(ix);
       GetCorrection(x,roc,dx);
       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
-      if (90.<=r0 && r0<=250.) {
+      if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
        Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
        h->SetBinContent(ix,iy,r1-r0);
       }
@@ -260,6 +305,7 @@ TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
        h->SetBinContent(ix,iy,0.);
     }
   }
+  delete tpcparam;
   return h;
 }
 
@@ -271,6 +317,8 @@ TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
   // The histogramm has nx times ny entries. 
   //
 
+  AliTPCParam* tpcparam = new AliTPCParamSR;
+
   TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
                     nx,-250.,250.,ny,-250.,250.);
   Float_t x[3],dx[3];
@@ -282,7 +330,7 @@ TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
       x[0]=h->GetXaxis()->GetBinCenter(ix);
       GetCorrection(x,roc,dx);
       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
-      if (90.<=r0 && r0<=250.) {
+      if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
        Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
        Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 
@@ -296,6 +344,39 @@ TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
        h->SetBinContent(ix,iy,0.);
     }
   }
+  delete tpcparam;
+  return h;
+}
+
+TH2F* AliTPCCorrection::CreateHistoDZinXY(Float_t z,Int_t nx,Int_t ny) {
+  //
+  // Simple plot functionality.
+  // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
+  // in respect to position z within the XY plane.
+  // The histogramm has nx times ny entries. 
+  //
+
+  AliTPCParam* tpcparam = new AliTPCParamSR;
+  TH2F *h=CreateTH2F("dz_xy",GetTitle(),"x [cm]","y [cm]","dz [cm]",
+                    nx,-250.,250.,ny,-250.,250.);
+  Float_t x[3],dx[3];
+  x[2]=z;
+  Int_t roc=z>0.?0:18; // FIXME
+  for (Int_t iy=1;iy<=ny;++iy) {
+    x[1]=h->GetYaxis()->GetBinCenter(iy);
+    for (Int_t ix=1;ix<=nx;++ix) {
+      x[0]=h->GetXaxis()->GetBinCenter(ix);
+      GetCorrection(x,roc,dx);
+      Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
+      if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
+       h->SetBinContent(ix,iy,dx[2]);
+      }
+      else
+       h->SetBinContent(ix,iy,0.);
+    }
+  }
+  delete tpcparam;
   return h;
 }
 
@@ -322,7 +403,6 @@ TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
       h->SetBinContent(iz,ir,r1-r0);
     }
   }
-  printf("SDF\n");
   return h;
 
 }
@@ -359,6 +439,32 @@ TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
   return h;
 }
 
+TH2F* AliTPCCorrection::CreateHistoDZinZR(Float_t phi,Int_t nz,Int_t nr) {
+  //
+  // Simple plot functionality.
+  // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz) 
+  // in respect to angle phi within the ZR plane.
+  // The histogramm has nx times ny entries. 
+  //
+  TH2F *h=CreateTH2F("dz_zr",GetTitle(),"z [cm]","r [cm]","dz [cm]",
+                    nz,-250.,250.,nr,85.,250.);
+  Float_t x[3],dx[3];
+  for (Int_t ir=1;ir<=nr;++ir) {
+    Float_t radius=h->GetYaxis()->GetBinCenter(ir);
+    x[0]=radius*TMath::Cos(phi);
+    x[1]=radius*TMath::Sin(phi);
+    for (Int_t iz=1;iz<=nz;++iz) {
+      x[2]=h->GetXaxis()->GetBinCenter(iz);
+      Int_t roc=x[2]>0.?0:18; // FIXME
+      GetCorrection(x,roc,dx);
+      h->SetBinContent(iz,ir,dx[2]);
+    }
+  }
+  return h;
+
+}
+
+
 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
                                   const char *xlabel,const char *ylabel,const char *zlabel,
                                  Int_t nbinsx,Double_t xlow,Double_t xup,
@@ -387,16 +493,14 @@ TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
   return h;
 }
 
-
 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
 
 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z, 
-                                                 const Double_t er[kNZ][kNR], Double_t &er_value )
-{
+                                                 const Double_t er[kNZ][kNR], Double_t &erValue ) {
   //
   // Interpolate table - 2D interpolation
   //
-  Double_t save_er[10] ;
+  Double_t saveEr[5] = {0,0,0,0,0};
 
   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
   Search( kNR,   fgkRList,  r,   fKLow   ) ;
@@ -406,16 +510,124 @@ void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double
   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 
   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
-      save_er[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
+      saveEr[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
   }
-  er_value = Interpolate( &fgkZList[fJLow], save_er, order, z )   ;
+  erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z )   ;
 
 }
 
+void AliTPCCorrection::Interpolate3DEdistortion( const Int_t order, const Double_t r, const Float_t phi, const Double_t z, 
+                                                const Double_t er[kNZ][kNPhi][kNR], const Double_t ephi[kNZ][kNPhi][kNR], const Double_t ez[kNZ][kNPhi][kNR],
+                                                Double_t &erValue, Double_t &ephiValue, Double_t &ezValue) {
+  //
+  // Interpolate table - 3D interpolation
+  //
+  
+  Double_t saveEr[5]= {0,0,0,0,0};
+  Double_t savedEr[5]= {0,0,0,0,0} ;
+
+  Double_t saveEphi[5]= {0,0,0,0,0};
+  Double_t savedEphi[5]= {0,0,0,0,0} ;
+
+  Double_t saveEz[5]= {0,0,0,0,0};
+  Double_t savedEz[5]= {0,0,0,0,0} ;
+
+  Search( kNZ,   fgkZList,   z,   fILow   ) ;
+  Search( kNPhi, fgkPhiList, z,   fJLow   ) ;
+  Search( kNR,   fgkRList,   r,   fKLow   ) ;
+
+  if ( fILow < 0 ) fILow = 0 ;   // check if out of range
+  if ( fJLow < 0 ) fJLow = 0 ;
+  if ( fKLow < 0 ) fKLow = 0 ;
+
+  if ( fILow + order  >=    kNZ - 1 ) fILow =   kNZ - 1 - order ;
+  if ( fJLow + order  >=  kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
+  if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
+
+  for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
+    for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
+      saveEr[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r )   ;
+      saveEphi[j-fJLow]   = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
+      saveEz[j-fJLow]     = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r )   ;
+    }
+    savedEr[i-fILow]     = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi )   ; 
+    savedEphi[i-fILow]   = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ; 
+    savedEz[i-fILow]     = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi )   ; 
+  }
+  erValue     = Interpolate( &fgkZList[fILow], savedEr, order, z )    ;
+  ephiValue   = Interpolate( &fgkZList[fILow], savedEphi, order, z )  ;
+  ezValue     = Interpolate( &fgkZList[fILow], savedEz, order, z )    ;
+
+}
+
+Double_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y, 
+                                             const Int_t nx,  const Int_t ny, const Double_t xv[], const Double_t yv[], 
+                                             const TMatrixD &array ) {
+  //
+  // Interpolate table (TMatrix format) - 2D interpolation
+  //
+
+  static  Int_t jlow = 0, klow = 0 ;
+  Double_t saveArray[5] = {0,0,0,0,0} ;
+
+  Search( nx,  xv,  x,   jlow  ) ;
+  Search( ny,  yv,  y,   klow  ) ;
+  if ( jlow < 0 ) jlow = 0 ;   // check if out of range
+  if ( klow < 0 ) klow = 0 ;
+  if ( jlow + order  >=    nx - 1 ) jlow =   nx - 1 - order ;
+  if ( klow + order  >=    ny - 1 ) klow =   ny - 1 - order ;
+
+  for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
+    {
+      Double_t *ajkl = &((TMatrixD&)array)(j,klow);
+      saveArray[j-jlow]  = Interpolate( &yv[klow], ajkl , order, y )   ;
+    }
+
+  return( Interpolate( &xv[jlow], saveArray, order, x ) )   ;
+
+}
+
+Double_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x,   const Double_t y,   const Double_t z,
+                                             const Int_t  nx,    const Int_t  ny,    const Int_t  nz,
+                                             const Double_t xv[], const Double_t yv[], const Double_t zv[],
+                                             TMatrixD **arrayofArrays ) {
+  //
+  // Interpolate table (TMatrix format) - 3D interpolation
+  //
+
+  static  Int_t ilow = 0, jlow = 0, klow = 0 ;
+  Double_t saveArray[5]= {0,0,0,0,0};
+  Double_t savedArray[5]= {0,0,0,0,0} ;
+
+  Search( nx, xv, x, ilow   ) ;
+  Search( ny, yv, y, jlow   ) ;
+  Search( nz, zv, z, klow   ) ;  
+
+  if ( ilow < 0 ) ilow = 0 ;   // check if out of range
+  if ( jlow < 0 ) jlow = 0 ;
+  if ( klow < 0 ) klow = 0 ;
+
+  if ( ilow + order  >=    nx - 1 ) ilow =   nx - 1 - order ;
+  if ( jlow + order  >=    ny - 1 ) jlow =   ny - 1 - order ;
+  if ( klow + order  >=    nz - 1 ) klow =   nz - 1 - order ;
+
+  for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
+    {
+      TMatrixD &table = *arrayofArrays[k] ;
+      for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
+       {
+         saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y )   ;
+       }
+      savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x )   ; 
+    }
+  return( Interpolate( &zv[klow], savedArray, order, z ) )   ;
+
+}
+
+
 
 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[], 
-                                      const Int_t order, const Double_t x )
-{
+                                      const Int_t order, const Double_t x ) {
   //
   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
   //
@@ -434,8 +646,7 @@ Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t
 }
 
 
-void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low )
-{
+void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
   //
   // Search an ordered table by starting at the most recently used point
   //
@@ -482,8 +693,590 @@ void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Dou
   
 }
 
+void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity, 
+                                          TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz, 
+                                          const Int_t rows, const Int_t columns, const Int_t iterations,
+                                          const Bool_t rocDisplacement ) {
+  //
+  // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
+  //
+  // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the 
+  // boundary conditions on the first and last rows, and the first and last columns.  The remainder of the 
+  // array can be blank or contain a preliminary guess at the solution.  The Charge density matrix contains 
+  // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if 
+  // you wish to solve Laplaces equation however it should not contain random numbers or you will get 
+  // random numbers back as a solution. 
+  // Poisson's equation is solved by iteratively relaxing the matrix to the final solution.  In order to 
+  // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution 
+  // space.  First it solves the problem on a very sparse grid by skipping rows and columns in the original 
+  // matrix.  Then it doubles the number of points and solves the problem again.  Then it doubles the 
+  // number of points and solves the problem again.  This happens several times until the maximum number
+  // of points has been included in the array.  
+  //
+  // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
+  // So rows == 2**M + 1 and columns == 2**N + 1.  The number of rows and columns can be different.
+  // 
+  // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
+  //
+  // Original code by Jim Thomas (STAR TPC Collaboration)
+  //
+
+  Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm; 
+
+  const Float_t  gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
+  const Float_t  gridSizeZ   =  fgkTPCZ0 / (columns-1) ;
+  const Float_t  ratio       =  gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
+
+  TMatrixD  arrayEr(rows,columns) ;
+  TMatrixD  arrayEz(rows,columns) ;
+
+  //Check that number of rows and columns is suitable for a binary expansion
+  
+  if ( !IsPowerOfTwo(rows-1) ) {
+    AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
+    return;
+  }
+  if ( !IsPowerOfTwo(columns-1) ) {
+    AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
+    return;
+  }
+  
+  // Solve Poisson's equation in cylindrical coordinates by relaxation technique
+  // Allow for different size grid spacing in R and Z directions
+  // Use a binary expansion of the size of the matrix to speed up the solution of the problem
+  
+  Int_t iOne = (rows-1)/4 ;
+  Int_t jOne = (columns-1)/4 ;
+  // Solve for N in 2**N, add one.
+  Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;  
+
+  for ( Int_t count = 0 ; count < loops ; count++ ) { 
+    // Loop while the matrix expands & the resolution increases.
+
+    Float_t tempGridSizeR = gridSizeR * iOne ;
+    Float_t tempRatio     = ratio * iOne * iOne / ( jOne * jOne ) ;
+    Float_t tempFourth    = 1.0 / (2.0 + 2.0*tempRatio) ;
+    
+    // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+    std::vector<float> coef1(rows) ;  
+    std::vector<float> coef2(rows) ;  
+
+    for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
+       Float_t radius = fgkIFCRadius + i*gridSizeR ;
+      coef1[i] = 1.0 + tempGridSizeR/(2*radius);
+      coef2[i] = 1.0 - tempGridSizeR/(2*radius);
+    }
+    
+    TMatrixD sumChargeDensity(rows,columns) ;
+
+    for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
+      Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
+      for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
+       if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
+       else {        
+         // Add up all enclosed charge density contributions within 1/2 unit in all directions
+         Float_t weight = 0.0 ;
+         Float_t sum    = 0.0 ;
+         sumChargeDensity(i,j) = 0.0 ;
+         for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
+           for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
+             if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
+             else
+               weight = 1.0 ;
+             // Note that this is cylindrical geometry
+             sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;  
+             sum += weight*radius ;
+           }
+         }
+         sumChargeDensity(i,j) /= sum ;
+       }
+        sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
+       }
+    }
+
+    for ( Int_t k = 1 ; k <= iterations; k++ ) {               
+      // Solve Poisson's Equation
+      // Over-relaxation index, must be >= 1 but < 2.  Arrange for it to evolve from 2 => 1 
+      // as interations increase.
+      Float_t overRelax   = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ; 
+      Float_t overRelaxM1 = overRelax - 1.0 ;
+      Float_t overRelaxtempFourth, overRelaxcoef5 ;
+      overRelaxtempFourth = overRelax * tempFourth ;
+      overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ; 
+
+      for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
+       for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
+
+         arrayV(i,j) = (   coef2[i]       *   arrayV(i-iOne,j)
+                         + tempRatio      * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
+                         - overRelaxcoef5 *   arrayV(i,j) 
+                         + coef1[i]       *   arrayV(i+iOne,j) 
+                         + sumChargeDensity(i,j) 
+                       ) * overRelaxtempFourth;
+       }
+      }
+
+      if ( k == iterations ) {    
+       // After full solution is achieved, copy low resolution solution into higher res array
+       for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
+         for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
+
+           if ( iOne > 1 ) {              
+             arrayV(i+iOne/2,j)                    =  ( arrayV(i+iOne,j) + arrayV(i,j)     ) / 2 ;
+             if ( i == iOne )  arrayV(i-iOne/2,j) =  ( arrayV(0,j)       + arrayV(iOne,j) ) / 2 ;
+           }
+           if ( jOne > 1 ) {
+             arrayV(i,j+jOne/2)                    =  ( arrayV(i,j+jOne) + arrayV(i,j) )     / 2 ;
+             if ( j == jOne )  arrayV(i,j-jOne/2) =  ( arrayV(i,0)       + arrayV(i,jOne) ) / 2 ;
+           }
+           if ( iOne > 1 && jOne > 1 ) {
+             arrayV(i+iOne/2,j+jOne/2) =  ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
+             if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
+             if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
+             // Note that this leaves a point at the upper left and lower right corners uninitialized. 
+             // -> Not a big deal.
+           }
+
+         }
+       }
+      }
+
+    }
+
+    iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
+    jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
+
+    sumChargeDensity.Clear();
+  }      
+
+  // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
+  for ( Int_t j = 0 ; j < columns ; j++ ) {      
+    for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
+    arrayEr(0,j)      =  -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;  
+    arrayEr(rows-1,j) =  -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ; 
+  }
+
+  // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
+  for ( Int_t i = 0 ; i < rows ; i++) {
+    for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
+    arrayEz(i,0)         =  -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;  
+    arrayEz(i,columns-1) =  -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ; 
+  }
+  
+  for ( Int_t i = 0 ; i < rows ; i++) {
+    // Note: go back and compare to old version of this code.  See notes below.
+    // JT Test ... attempt to divide by real Ez not Ez to first order
+    for ( Int_t j = 0 ; j < columns ; j++ ) {
+      arrayEz(i,j) += ezField;
+      // This adds back the overall Z gradient of the field (main E field component)
+    } 
+    // Warning: (-=) assumes you are using an error potetial without the overall Field included
+  }                                 
+  
+  // Integrate Er/Ez from Z to zero
+  for ( Int_t j = 0 ; j < columns ; j++ )  {     
+    for ( Int_t i = 0 ; i < rows ; i++ ) {
+      
+      Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.  
+      arrayErOverEz(i,j) = 0.0 ;
+      arrayDeltaEz(i,j) = 0.0 ;
+      
+      for ( Int_t k = j ; k < columns ; k++ ) {
+       arrayErOverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
+       arrayDeltaEz(i,j)   +=  index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
+       if ( index != 4 )  index = 4; else index = 2 ;
+      }
+      if ( index == 4 ) {
+       arrayErOverEz(i,j)  -=  (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
+       arrayDeltaEz(i,j)   -=  (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
+      }
+      if ( index == 2 ) {
+       arrayErOverEz(i,j)  +=  (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2) 
+                                                   -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
+       arrayDeltaEz(i,j)   +=  (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField) 
+                                                   -2.5*(arrayEz(i,columns-1)-ezField));
+      }
+      if ( j == columns-2 ) {
+       arrayErOverEz(i,j) =  (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
+                                                 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
+       arrayDeltaEz(i,j)  =  (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
+                                                 +1.5*(arrayEz(i,columns-1)-ezField) ) ;
+      }
+      if ( j == columns-1 ) {
+       arrayErOverEz(i,j) =  0.0 ;
+       arrayDeltaEz(i,j)  =  0.0 ;
+      }
+    }
+  }
+  
+  // calculate z distortion from the integrated Delta Ez residuals
+  // and include the aquivalence (Volt to cm) of the ROC shift !!
+
+  for ( Int_t j = 0 ; j < columns ; j++ )  {     
+    for ( Int_t i = 0 ; i < rows ; i++ ) {
+
+      // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
+      arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
+
+      // ROC Potential in cm aquivalent
+      Double_t dzROCShift =  arrayV(i, columns -1)/ezField;  
+      if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift;  // add the ROC misaligment
+
+    }
+  }
+  arrayEr.Clear();
+  arrayEz.Clear();
+
+}
+
+void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities, 
+                   TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
+                   const Int_t rows, const Int_t columns,  const Int_t phislices, 
+                   const Float_t deltaphi, const Int_t iterations, const Int_t symmetry,
+                   Bool_t rocDisplacement  ) {
+  //
+  // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
+  //
+  //    NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.  
+  //    The number of rows and COLUMNS can be different.
+  //
+  //    ROWS       ==  2**M + 1  
+  //    COLUMNS    ==  2**N + 1  
+  //    PHISLICES  ==  Arbitrary but greater than 3
+  //
+  //    DeltaPhi in Radians
+  //
+  //    SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
+  //             = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
+  //
+  // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
+
+  const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm; 
+
+  const Float_t  gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
+  const Float_t  gridSizePhi =  deltaphi ;
+  const Float_t  gridSizeZ   =  fgkTPCZ0 / (columns-1) ;
+  const Float_t  ratioPhi    =  gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
+  const Float_t  ratioZ      =  gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
+
+  TMatrixD arrayE(rows,columns) ;
+
+  // Check that the number of rows and columns is suitable for a binary expansion
+  if ( !IsPowerOfTwo((rows-1))    ) {  
+    AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1"); 
+    return; }
+  if ( !IsPowerOfTwo((columns-1)) ) { 
+    AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
+    return; }
+  if ( phislices <= 3   )  { 
+    AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
+    return; }
+  if  ( phislices > 1000 ) { 
+    AliError("Poisson3D  phislices > 1000 is not allowed (nor wise) ");  
+    return; }  
+  
+  // Solve Poisson's equation in cylindrical coordinates by relaxation technique
+  // Allow for different size grid spacing in R and Z directions
+  // Use a binary expansion of the matrix to speed up the solution of the problem
+
+  Int_t loops, mplus, mminus, signplus, signminus  ;
+  Int_t ione = (rows-1)/4 ;
+  Int_t jone = (columns-1)/4 ;
+  loops = TMath::Max(ione, jone) ;      // Calculate the number of loops for the binary expansion
+  loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ;  // Solve for N in 2**N
+
+  TMatrixD* arrayofSumChargeDensities[1000] ;    // Create temporary arrays to store low resolution charge arrays
+
+  for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
+
+  for ( Int_t count = 0 ; count < loops ; count++ ) {      // START the master loop and do the binary expansion
+   
+    Float_t  tempgridSizeR   =  gridSizeR  * ione ;
+    Float_t  tempratioPhi    =  ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
+    Float_t  tempratioZ      =  ratioZ   * ione * ione / ( jone * jone ) ;
+
+    std::vector<float> coef1(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+    std::vector<float> coef2(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+    std::vector<float> coef3(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+    std::vector<float> coef4(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+
+    for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
+      Float_t radius = fgkIFCRadius + i*gridSizeR ;
+      coef1[i] = 1.0 + tempgridSizeR/(2*radius);
+      coef2[i] = 1.0 - tempgridSizeR/(2*radius);
+      coef3[i] = tempratioPhi/(radius*radius);
+      coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
+    }
+
+    for ( Int_t m = 0 ; m < phislices ; m++ ) {
+      TMatrixD &chargeDensity    = *arrayofChargeDensities[m] ;
+      TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
+      for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
+       Float_t radius = fgkIFCRadius + i*gridSizeR ;
+       for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
+         if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
+         else {           // Add up all enclosed charge density contributions within 1/2 unit in all directions
+           Float_t weight = 0.0 ;
+           Float_t sum    = 0.0 ;
+           sumChargeDensity(i,j) = 0.0 ;
+           for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
+             for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
+               if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
+               else
+                 weight = 1.0 ; 
+               sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;  
+               sum += weight*radius ;
+             }
+           }
+           sumChargeDensity(i,j) /= sum ;
+         }
+          sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
+       }
+      }
+    }
+
+    for ( Int_t k = 1 ; k <= iterations; k++ ) {
+
+      // over-relaxation index, >= 1 but < 2
+      Float_t overRelax   = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ; 
+      Float_t overRelaxM1 = overRelax - 1.0 ;
+
+      std::vector<float> overRelaxcoef4(rows) ;  // Do this the standard C++ way to avoid gcc extensions
+      std::vector<float> overRelaxcoef5(rows) ;  // Do this the standard C++ way to avoid gcc extensions
+
+      for ( Int_t i = ione ; i < rows-1 ; i+=ione ) { 
+       overRelaxcoef4[i] = overRelax * coef4[i] ;
+       overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ; 
+      }
+
+      for ( Int_t m = 0 ; m < phislices ; m++ ) {
+
+       mplus  = m + 1;   signplus  = 1 ; 
+       mminus = m - 1 ;  signminus = 1 ;
+       if (symmetry==1) {  // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
+         if ( mplus  > phislices-1 ) mplus  = phislices - 2 ;
+         if ( mminus < 0 )           mminus = 1 ;
+       }
+       else if (symmetry==-1) {   // Anti-symmetry in phi
+         if ( mplus  > phislices-1 ) { mplus  = phislices - 2 ; signplus  = -1 ; }
+         if ( mminus < 0 )           { mminus = 1 ;             signminus = -1 ; } 
+       }
+               else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
+         if ( mplus  > phislices-1 ) mplus  = m + 1 - phislices ;
+         if ( mminus < 0 )           mminus = m - 1 + phislices ;
+       }
+       TMatrixD& arrayV    =  *arrayofArrayV[m] ;
+       TMatrixD& arrayVP   =  *arrayofArrayV[mplus] ;
+       TMatrixD& arrayVM   =  *arrayofArrayV[mminus] ;
+       TMatrixD& sumChargeDensity =  *arrayofSumChargeDensities[m] ;
+
+       for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
+         for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
+
+            arrayV(i,j) = (   coef2[i]          *   arrayV(i-ione,j)
+                           + tempratioZ        * ( arrayV(i,j-jone)  +  arrayV(i,j+jone) )
+                           - overRelaxcoef5[i] *   arrayV(i,j) 
+                           + coef1[i]          *   arrayV(i+ione,j)  
+                           + coef3[i]          * ( signplus*arrayVP(i,j)       +  signminus*arrayVM(i,j) )
+                           + sumChargeDensity(i,j) 
+                         ) * overRelaxcoef4[i] ;     
+           // Note: over-relax the solution at each step.  This speeds up the convergance.
+
+         }
+       }
+
+       if ( k == iterations ) {   // After full solution is achieved, copy low resolution solution into higher res array
+         for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
+           for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
+             
+             if ( ione > 1 ) {              
+               arrayV(i+ione/2,j)                    =  ( arrayV(i+ione,j) + arrayV(i,j)     ) / 2 ;
+               if ( i == ione )  arrayV(i-ione/2,j) =  ( arrayV(0,j)       + arrayV(ione,j) ) / 2 ;
+             }
+             if ( jone > 1 ) {
+               arrayV(i,j+jone/2)                    =  ( arrayV(i,j+jone) + arrayV(i,j) )     / 2 ;
+               if ( j == jone )  arrayV(i,j-jone/2) =  ( arrayV(i,0)       + arrayV(i,jone) ) / 2 ;
+             }
+             if ( ione > 1 && jone > 1 ) {
+               arrayV(i+ione/2,j+jone/2) =  ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
+               if ( i == ione ) arrayV(i-ione/2,j-jone/2) =   ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
+               if ( j == jone ) arrayV(i-ione/2,j-jone/2) =   ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
+               // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
+             }
+           }       
+         }
+       }
+
+      }
+    }      
+
+    ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
+    jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
+
+  }
+  
+  //Differentiate V(r) and solve for E(r) using special equations for the first and last row
+  //Integrate E(r)/E(z) from point of origin to pad plane
+
+  for ( Int_t m = 0 ; m < phislices ; m++ ) {
+    TMatrixD& arrayV    =  *arrayofArrayV[m] ;
+    TMatrixD& eroverEz  =  *arrayofEroverEz[m] ;
+    
+    for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {  // Count backwards to facilitate integration over Z
+      
+      // Differentiate in R
+      for ( Int_t i = 1 ; i < rows-1 ; i++ )  arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
+      arrayE(0,j)      =  -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;  
+      arrayE(rows-1,j) =  -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ; 
+      // Integrate over Z
+      for ( Int_t i = 0 ; i < rows ; i++ ) {
+       Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.  
+       eroverEz(i,j) = 0.0 ;
+       for ( Int_t k = j ; k < columns ; k++ ) {
+         
+         eroverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
+         if ( index != 4 )  index = 4; else index = 2 ;
+       }
+       if ( index == 4 ) eroverEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
+       if ( index == 2 ) eroverEz(i,j)  +=  
+         (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
+       if ( j == columns-2 ) eroverEz(i,j) =  
+         (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
+       if ( j == columns-1 ) eroverEz(i,j) =  0.0 ;
+      }
+    }
+    // if ( m == 0 ) { TCanvas*  c1 =  new TCanvas("erOverEz","erOverEz",50,50,840,600) ;  c1 -> cd() ;
+    // eroverEz.Draw("surf") ; } // JT test
+  }
+  
+  //Differentiate V(r) and solve for E(phi) 
+  //Integrate E(phi)/E(z) from point of origin to pad plane
+
+  for ( Int_t m = 0 ; m < phislices ; m++ ) {
+    
+    mplus  = m + 1;   signplus  = 1 ; 
+    mminus = m - 1 ;  signminus = 1 ; 
+    if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
+      if ( mplus  > phislices-1 ) mplus  = phislices - 2 ;
+      if ( mminus < 0 )           mminus = 1 ;
+    }
+    else if (symmetry==-1) {       // Anti-symmetry in phi
+      if ( mplus  > phislices-1 ) { mplus  = phislices - 2 ;  signplus  = -1 ; }
+      if ( mminus < 0 )           { mminus = 1 ;                   signminus = -1 ; } 
+    }
+    else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
+      if ( mplus  > phislices-1 ) mplus  = m + 1 - phislices ;
+      if ( mminus < 0 )           mminus = m - 1 + phislices ;
+    }
+    TMatrixD &arrayVP     =  *arrayofArrayV[mplus] ;
+    TMatrixD &arrayVM     =  *arrayofArrayV[mminus] ;
+    TMatrixD &ePhioverEz  =  *arrayofEPhioverEz[m] ;
+    for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
+      // Differentiate in Phi
+      for ( Int_t i = 0 ; i < rows ; i++ ) {
+       Float_t radius = fgkIFCRadius + i*gridSizeR ;
+       arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
+      }
+      // Integrate over Z
+      for ( Int_t i = 0 ; i < rows ; i++ ) {
+       Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.  
+       ePhioverEz(i,j) = 0.0 ;
+       for ( Int_t k = j ; k < columns ; k++ ) {
+         
+         ePhioverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
+         if ( index != 4 )  index = 4; else index = 2 ;
+       }
+       if ( index == 4 ) ePhioverEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
+       if ( index == 2 ) ePhioverEz(i,j)  +=  
+         (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
+       if ( j == columns-2 ) ePhioverEz(i,j) =  
+         (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
+       if ( j == columns-1 ) ePhioverEz(i,j) =  0.0 ;
+      }
+    }
+    // if ( m == 5 ) { TCanvas* c2 =  new TCanvas("arrayE","arrayE",50,50,840,600) ;  c2 -> cd() ;
+    // arrayE.Draw("surf") ; } // JT test
+  }
+  
+
+  // Differentiate V(r) and solve for E(z) using special equations for the first and last row
+  // Integrate (E(z)-Ezstd) from point of origin to pad plane
+  
+  for ( Int_t m = 0 ; m < phislices ; m++ ) {
+    TMatrixD& arrayV   =  *arrayofArrayV[m] ;
+    TMatrixD& deltaEz  =  *arrayofDeltaEz[m] ;
+    
+    // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
+    for ( Int_t i = 0 ; i < rows ; i++) {
+      for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
+      arrayE(i,0)         =  -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;  
+      arrayE(i,columns-1) =  -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ; 
+    }
+    
+    for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {  // Count backwards to facilitate integration over Z
+      // Integrate over Z
+      for ( Int_t i = 0 ; i < rows ; i++ ) {
+       Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.  
+       deltaEz(i,j) = 0.0 ;
+       for ( Int_t k = j ; k < columns ; k++ ) {
+         deltaEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k) ;
+         if ( index != 4 )  index = 4; else index = 2 ;
+       }
+       if ( index == 4 ) deltaEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1) ;
+       if ( index == 2 ) deltaEz(i,j)  +=  
+         (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
+       if ( j == columns-2 ) deltaEz(i,j) =  
+         (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
+       if ( j == columns-1 ) deltaEz(i,j) =  0.0 ;
+      }
+    }
+    // if ( m == 0 ) { TCanvas*  c1 =  new TCanvas("erOverEz","erOverEz",50,50,840,600) ;  c1 -> cd() ;
+    // eroverEz.Draw("surf") ; } // JT test
+    
+    // calculate z distortion from the integrated Delta Ez residuals
+    // and include the aquivalence (Volt to cm) of the ROC shift !!
+    
+    for ( Int_t j = 0 ; j < columns ; j++ )  {   
+      for ( Int_t i = 0 ; i < rows ; i++ ) {
+       
+       // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
+       deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
+       
+       // ROC Potential in cm aquivalent
+       Double_t dzROCShift =  arrayV(i, columns -1)/ezField;  
+       if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift;  // add the ROC misaligment
+       
+      }
+    }
+
+  } // end loop over phi
+  
 
-AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir,TTreeSRedirector *pcstream){
+  for ( Int_t k = 0 ; k < phislices ; k++ )
+    {
+      arrayofSumChargeDensities[k]->Delete() ;
+    }
+  
+
+
+  arrayE.Clear();
+}
+
+
+Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
+  //
+  // Helperfunction: Check if integer is a power of 2
+  //
+  Int_t j = 0;
+  while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
+  if ( j == 1 ) return(1) ;  // True
+  return(0) ;                // False
+}
+
+
+AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
   //
   // Fit the track parameters - without and with distortion
   // 1. Space points in the TPC are simulated along the trajectory  
@@ -502,31 +1295,38 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara
   // track1.fP[2] - sinus of local inclination angle
   // track1.fP[3] - tangent of deep angle
   // track1.fP[4] - 1/pt
+
   AliTPCROC * roc = AliTPCROC::Instance();
   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
-    
   const Double_t kMaxSnp = 0.85;  
   const Double_t kSigmaY=0.1;
   const Double_t kSigmaZ=0.1;
+  const Double_t kMaxR=500;
+  const Double_t kMaxZ=500;
   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+  Int_t npoints1=0;
+  Int_t npoints2=0;
 
   AliExternalTrackParam  track(trackIn); // 
   // generate points
   AliTrackPointArray pointArray0(npoints0);
   AliTrackPointArray pointArray1(npoints0);
   Double_t xyz[3];
-  AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
+  if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
   //
   // simulate the track
   Int_t npoints=0;
   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
-    AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp);
+    if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
     track.GetXYZ(xyz);
-    xyz[0]+=gRandom->Gaus(0,0.001);
-    xyz[1]+=gRandom->Gaus(0,0.001);
+    xyz[0]+=gRandom->Gaus(0,0.00005);
+    xyz[1]+=gRandom->Gaus(0,0.00005);
+    xyz[2]+=gRandom->Gaus(0,0.00005);
+    if (TMath::Abs(track.GetZ())>kMaxZ) break;
+    if (TMath::Abs(track.GetX())>kMaxR) break;
     AliTrackPoint pIn0;                               // space point          
     AliTrackPoint pIn1;
     Int_t sector= (xyz[2]>0)? 0:18;
@@ -550,6 +1350,7 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara
     npoints++;
     if (npoints>=npoints0) break;
   }
+  if (npoints<npoints0/2) return 0;
   //
   // refit track
   //
@@ -569,7 +1370,6 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara
   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
 
-
   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
     Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
     //
@@ -580,8 +1380,13 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara
     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
     //
-    AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
-    AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
+    if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
+    if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
+    if (TMath::Abs(track0->GetZ())>kMaxZ) break;
+    if (TMath::Abs(track0->GetX())>kMaxR) break;
+    if (TMath::Abs(track1->GetZ())>kMaxZ) break;
+    if (TMath::Abs(track1->GetX())>kMaxR) break;
+
     track.GetXYZ(xyz);  // distorted track also propagated to the same reference radius
     //
     Double_t pointPos[2]={0,0};
@@ -591,23 +1396,25 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara
     pointCov[0]=prot0.GetCov()[3];//simay^2
     pointCov[1]=prot0.GetCov()[4];//sigmayz
     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
-    track0->Update(pointPos,pointCov);
+    if (!track0->Update(pointPos,pointCov)) break;
     //
     Double_t deltaX=prot1.GetX()-prot0.GetX();   // delta X 
     Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp()));  // deltaY due  delta X
     Double_t deltaZX=deltaX*track1->GetTgl();                           // deltaZ due  delta X
 
-    pointPos[0]=prot1.GetY()-deltaYX;//local y
-    pointPos[1]=prot1.GetZ()-deltaZX;//local z
+    pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
+    pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
     pointCov[0]=prot1.GetCov()[3];//simay^2
     pointCov[1]=prot1.GetCov()[4];//sigmayz
     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
-    track1->Update(pointPos,pointCov);
+    if (!track1->Update(pointPos,pointCov)) break;
+    npoints1++;
+    npoints2++;
   }
-
+  if (npoints2<npoints)  return 0;
   AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
   track1->Rotate(track0->GetAlpha());
-  track1->PropagateTo(track0->GetX(),AliTrackerBase::GetBz());
+  AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
 
   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
     "point0.="<<&pointArray0<<   //  points
@@ -688,7 +1495,7 @@ TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
 
 
 
-void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, TObjArray * corrArray, Int_t step, Bool_t debug ){
+void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
   //
   // Make a fit tree:
   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
@@ -705,10 +1512,11 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
   // corrArray - array with partial corrections
   // step      - skipe entries  - if 1 all entries processed - it is slow
   // debug     0 if debug on also space points dumped - it is slow
+
   const Double_t kMaxSnp = 0.85;  
   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
   //  const Double_t kB2C=-0.299792458e-3;
-  const Int_t kMinEntries=50;
+  const Int_t kMinEntries=50; 
   Double_t phi,theta, snp, mean,rms, entries;
   tinput->SetBranchAddress("theta",&theta);
   tinput->SetBranchAddress("phi", &phi);
@@ -720,17 +1528,19 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
   //
   Int_t nentries=tinput->GetEntries();
   Int_t ncorr=corrArray->GetEntries();
-  Double_t corrections[100]; //
+  Double_t corrections[100]={0}; //
   Double_t tPar[5];
   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
   Double_t refX=0;
   Int_t dir=0;
-  if (dtype==0) {refX=85; dir=-1;}
-  if (dtype==1) {refX=275; dir=1;}
-  if (dtype==2) {refX=85; dir=-1;}
+  if (dtype==0) {refX=85.; dir=-1;}
+  if (dtype==1) {refX=275.; dir=1;}
+  if (dtype==2) {refX=85.; dir=-1;}
+  if (dtype==3) {refX=360.; dir=-1;}
   //
   for (Int_t ientry=0; ientry<nentries; ientry+=step){
     tinput->GetEntry(ientry);
+    if (TMath::Abs(snp)>kMaxSnp) continue;
     tPar[0]=0;
     tPar[1]=theta*refX;
     tPar[2]=snp;
@@ -739,7 +1549,12 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
     Double_t bz=AliTrackerBase::GetBz();
     if (refX>10. && TMath::Abs(bz)>0.1 )  tPar[4]=snp/(refX*bz*kB2C*2);
     tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
-    if (TMath::Abs(snp)>0.251) continue;
+    AliExternalTrackParam track(refX,phi,tPar,cov);
+    Double_t xyz[3];
+    track.GetXYZ(xyz);
+    Int_t id=0;
+    Double_t dRrec=0; // dummy value - needed for points - e.g for laser
+    if (ptype==4 &&bz<0) mean*=-1;  // interpret as curvature
     (*pcstream)<<"fit"<<
       "bz="<<bz<<         // magnetic filed used
       "dtype="<<dtype<<   // detector match type
@@ -749,34 +1564,43 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
       "snp="<<snp<<       // snp
       "mean="<<mean<<     // mean dist value
       "rms="<<rms<<       // rms
+      "gx="<<xyz[0]<<         // global position at reference
+      "gy="<<xyz[1]<<         // global position at reference
+      "gz="<<xyz[2]<<         // global position at reference  
+      "dRrec="<<dRrec<<      // delta Radius in reconstruction
+      "id="<<id<<             // track id
       "entries="<<entries;// number of entries in bin
     //
     for (Int_t icorr=0; icorr<ncorr; icorr++) {
       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
       corrections[icorr]=0;
       if (entries>kMinEntries){
-       if (dtype==0) {
-         refX=85; dir=-1;
-       }
-       if (dtype==1) {
-         refX=275; dir=1;
-       }
-       if (dtype==2) {
-         refX=0; dir=-1;
-       }
-       //
        AliExternalTrackParam trackIn(refX,phi,tPar,cov);
        AliExternalTrackParam *trackOut = 0;
        if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
        if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
-       AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kFALSE,kMaxSnp);
-       AliTrackerBase::PropagateTrackToBxByBz(trackOut,refX,kMass,3,kFALSE,kMaxSnp);
+       if (dtype==0) {refX=85.; dir=-1;}
+       if (dtype==1) {refX=275.; dir=1;}
+       if (dtype==2) {refX=0; dir=-1;}
+       if (dtype==3) {refX=360.; dir=-1;}
        //
-       corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
-       delete trackOut;      
+       if (trackOut){
+         AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
+         trackOut->Rotate(trackIn.GetAlpha());
+         trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
+         //
+         corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
+         delete trackOut;      
+       }else{
+         corrections[icorr]=0;
+       }
+       if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature
       }      
+      Double_t dRdummy=0;
       (*pcstream)<<"fit"<<
-       Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
+       Form("%s=",corr->GetName())<<corrections[icorr]<<   // dump correction value
+       Form("dR%s=",corr->GetName())<<dRdummy;   // dump dummy correction value not needed for tracks 
+                                                  // for points it is neccessary
     }
     (*pcstream)<<"fit"<<"\n";
   }
@@ -785,8 +1609,130 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
 
 
 
+void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
+  //
+  // Make a laser fit tree for global minimization
+  //
+  const Double_t cutErrY=0.1;
+  const Double_t cutErrZ=0.1;
+  const Double_t kEpsilon=0.00000001;
+  TVectorD *vecdY=0;
+  TVectorD *vecdZ=0;
+  TVectorD *veceY=0;
+  TVectorD *veceZ=0;
+  AliTPCLaserTrack *ltr=0;
+  AliTPCLaserTrack::LoadTracks();
+  tree->SetBranchAddress("dY.",&vecdY);
+  tree->SetBranchAddress("dZ.",&vecdZ);
+  tree->SetBranchAddress("eY.",&veceY);
+  tree->SetBranchAddress("eZ.",&veceZ);
+  tree->SetBranchAddress("LTr.",&ltr);
+  Int_t entries= tree->GetEntries();
+  TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
+  Double_t bz=AliTrackerBase::GetBz();
+  // 
+
+  for (Int_t ientry=0; ientry<entries; ientry++){
+    tree->GetEntry(ientry);
+    if (!ltr->GetVecGX()){
+      ltr->UpdatePoints();
+    }
+    TVectorD * delta= (itype==0)? vecdY:vecdZ;
+    TVectorD * err= (itype==0)? veceY:veceZ;
+    
+    for (Int_t irow=0; irow<159; irow++){
+      Int_t nentries = 1000;
+      if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
+      if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
+      Int_t dtype=4;
+      Double_t phi   =(*ltr->GetVecPhi())[irow];
+      Double_t theta =ltr->GetTgl();
+      Double_t mean=delta->GetMatrixArray()[irow];
+      Double_t gx=0,gy=0,gz=0;
+      Double_t snp = (*ltr->GetVecP2())[irow];
+      Double_t rms = 0.1+err->GetMatrixArray()[irow];
+      gx = (*ltr->GetVecGX())[irow];
+      gy = (*ltr->GetVecGY())[irow];
+      gz = (*ltr->GetVecGZ())[irow];
+      Int_t bundle= ltr->GetBundle();
+      Double_t dRrec=0;
+      //
+      // get delta R used in reconstruction
+      AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();  
+      AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
+      const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
+      Double_t xyz0[3]={gx,gy,gz};
+      Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
+      //
+      // old ExB correction 
+      //      
+      if(recoParam&&recoParam->GetUseExBCorrection()) {        
+       Double_t xyz1[3]={gx,gy,gz};
+       calib->GetExB()->Correct(xyz0,xyz1);
+       Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
+       dRrec=oldR-newR;
+      } 
+      if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
+       Float_t xyz1[3]={gx,gy,gz};
+       Int_t sector=(gz>0)?0:18;
+       correction->CorrectPoint(xyz1, sector);
+       Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
+       dRrec=oldR-newR;
+      } 
+
+
+      (*pcstream)<<"fit"<<
+       "bz="<<bz<<         // magnetic filed used
+       "dtype="<<dtype<<   // detector match type
+       "ptype="<<itype<<   // parameter type
+       "theta="<<theta<<   // theta
+       "phi="<<phi<<       // phi 
+       "snp="<<snp<<       // snp
+       "mean="<<mean<<     // mean dist value
+       "rms="<<rms<<       // rms
+       "gx="<<gx<<         // global position
+       "gy="<<gy<<         // global position
+       "gz="<<gz<<         // global position
+       "dRrec="<<dRrec<<      // delta Radius in reconstruction
+       "id="<<bundle<<     //bundle
+       "entries="<<nentries;// number of entries in bin
+      //
+      //    
+      Double_t ky = TMath::Tan(TMath::ASin(snp));
+      Int_t ncorr = corrArray->GetEntries();
+      Double_t r0   = TMath::Sqrt(gx*gx+gy*gy);
+      Double_t phi0 = TMath::ATan2(gy,gx);
+      Double_t distortions[1000]={0};
+      Double_t distortionsR[1000]={0};
+      for (Int_t icorr=0; icorr<ncorr; icorr++) {
+       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
+       Float_t distPoint[3]={gx,gy,gz}; 
+       Int_t sector= (gz>0)? 0:18;
+       if (r0>80){
+         corr->DistortPoint(distPoint, sector);
+       }
+       // Double_t value=distPoint[2]-gz;
+       if (itype==0){
+         Double_t r1   = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
+         Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
+         Double_t drphi= r0*(phi1-phi0);
+         Double_t dr   = r1-r0;
+         distortions[icorr]  = drphi-ky*dr;
+         distortionsR[icorr] = dr;
+       }
+       (*pcstream)<<"fit"<<
+         Form("%s=",corr->GetName())<<distortions[icorr]<<    // dump correction value
+         Form("dR%s=",corr->GetName())<<distortionsR[icorr];   // dump correction R  value
+      }
+      (*pcstream)<<"fit"<<"\n";
+    }
+  }
+  delete pcstream;
+}
+
+
 
-void   AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector *pcstream, const char* hname, Int_t run){
+void   AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
   //
   // make a distortion map out ou fthe residual histogram
   // Results are written to the debug streamer - pcstream
@@ -882,3 +1828,158 @@ void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *com
   gStorage->Put(this, (*id1), metaData);
 }
 
+
+void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
+  //
+  // Fast method to simulate the influence of the given distortion on the vertex reconstruction
+  //
+
+  AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
+  if (!magF) AliError("Magneticd field - not initialized");
+  Double_t bz = magF->SolenoidField(); //field in kGauss
+  printf("bz: %lf\n",bz);
+  AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
+
+  TObjArray   aTrk;              // Original Track array of Aside
+  TObjArray   daTrk;             // Distorted Track array of A side
+  UShort_t    *aId = new UShort_t[nTracks];      // A side Track ID
+  TObjArray   cTrk;               
+  TObjArray   dcTrk;
+  UShort_t    *cId = new UShort_t [nTracks];
+  Int_t id=0; 
+  Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+  TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
+  fpt.SetParameters(7.24,0.120);
+  fpt.SetNpx(10000);
+  for(Int_t nt=0; nt<nTracks; nt++){
+    Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
+    Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
+    Double_t pt = fpt.GetRandom(); // momentum for f1
+    //   printf("phi %lf  eta %lf pt %lf\n",phi,eta,pt);
+    Short_t sign=1;
+    if(gRandom->Rndm() < 0.5){
+      sign =1;
+    }else{
+      sign=-1;
+    }
+
+    Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
+    Double_t pxyz[3];
+    pxyz[0]=pt*TMath::Cos(phi);
+    pxyz[1]=pt*TMath::Sin(phi);
+    pxyz[2]=pt*TMath::Tan(theta);
+    Double_t cv[21]={0};
+    AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
+
+    Double_t refX=1.;
+    Int_t dir=-1;
+    AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir,  NULL);
+    if (!td) continue;
+    if (pcstream) (*pcstream)<<"track"<<
+      "eta="<<eta<<
+      "theta="<<theta<<
+      "tOrig.="<<t<<
+      "td.="<<td<<
+      "\n";
+    if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
+      if (td){
+       daTrk.AddLast(td);
+       aTrk.AddLast(t);
+       Int_t nn=aTrk.GetEntriesFast();
+       aId[nn]=id;
+      }
+    }else if(( eta<-0.07 )&&( eta>-etaCuts )){
+      if (td){
+       dcTrk.AddLast(td);
+       cTrk.AddLast(t);
+       Int_t nn=cTrk.GetEntriesFast();
+       cId[nn]=id;
+      }
+    }
+    id++;  
+  }// end of track loop
+
+  vertexer->SetTPCMode();
+  vertexer->SetConstraintOff();
+
+  aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));  
+  avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
+  cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));  
+  cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
+  if (pcstream) (*pcstream)<<"vertex"<<
+    "x="<<orgVertex[0]<<
+    "y="<<orgVertex[1]<<
+    "z="<<orgVertex[2]<<
+    "av.="<<&aV<<              // distorted vertex A side
+    "cv.="<<&cV<<              // distroted vertex C side
+    "avO.="<<&avOrg<<         // original vertex A side
+    "cvO.="<<&cvOrg<<
+    "\n";
+  delete []aId;
+  delete []cId;
+}
+
+void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
+  //
+  // make correction available for visualization using 
+  // TFormula, TFX and TTree::Draw 
+  // important in order to check corrections and also compute dervied variables 
+  // e.g correction partial derivatives
+  //
+  // NOTE - class is not owner of correction
+  //     
+  if (!fgVisualCorrection) fgVisualCorrection=new TObjArray;
+  if (position!=0&&position>=fgVisualCorrection->GetEntriesFast())
+    fgVisualCorrection->Expand(position*2);
+  fgVisualCorrection->AddAt(corr, position);
+}
+
+
+
+Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t r, Double_t kZ, Int_t axisType, Int_t corrType){
+  //
+  // calculate the correction at given position - check the geffCorr
+  //
+  if (!fgVisualCorrection) return 0;
+  AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
+  if (!corr) return 0;
+
+  Double_t phi=sector*TMath::Pi()/9.;
+  Double_t gx = r*TMath::Cos(phi);
+  Double_t gy = r*TMath::Sin(phi);
+  Double_t gz = r*kZ;
+  Int_t nsector=(gz>0) ? 0:18; 
+  //
+  //
+  //
+  Float_t distPoint[3]={gx,gy,gz};
+  corr->DistortPoint(distPoint, nsector);
+  Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
+  Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
+  Double_t phi0=TMath::ATan2(gy,gx);
+  Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
+  if (axisType==0) return r1-r0;
+  if (axisType==1) return (phi1-phi0)*r0;
+  if (axisType==2) return distPoint[2]-gz;
+  return phi1-phi0;
+}
+
+Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
+  //
+  // return correction at given x,y,z
+  // 
+  if (!fgVisualCorrection) return 0;
+  AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
+  if (!corr) return 0;
+  Double_t phi0= TMath::ATan2(gy,gx);
+  Int_t nsector=(gz>0) ? 0:18; 
+  Float_t distPoint[3]={gx,gy,gz};
+  corr->DistortPoint(distPoint, nsector);
+  Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
+  Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
+  Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
+  if (axisType==0) return r1-r0;
+  if (axisType==1) return (phi1-phi0)*r0;
+  if (axisType==2) return distPoint[2]-gz;
+  return phi1-phi0;
+}