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