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[10] ;
 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[10],   savedEr[10] ;
 527   Double_t saveEphi[10], savedEphi[10] ;
 528   Double_t saveEz[10],   savedEz[10] ;
 529
 530   Search( kNZ,   fgkZList,   z,   fILow   ) ;
 531   Search( kNPhi, fgkPhiList, z,   fJLow   ) ;
 532   Search( kNR,   fgkRList,   r,   fKLow   ) ;
 533
 534   if ( fILow < 0 ) fILow = 0 ;   // check if out of range
 535   if ( fJLow < 0 ) fJLow = 0 ;
 536   if ( fKLow < 0 ) fKLow = 0 ;
 537
 538   if ( fILow + order  >=    kNZ - 1 ) fILow =   kNZ - 1 - order ;
 539   if ( fJLow + order  >=  kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
 540   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 541
 542   for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
 543     for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
 544       saveEr[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r )   ;
 545       saveEphi[j-fJLow]   = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
 546       saveEz[j-fJLow]     = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r )   ;
 547     }
 548     savedEr[i-fILow]     = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi )   ;
 549     savedEphi[i-fILow]   = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
 550     savedEz[i-fILow]     = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi )   ;
 551   }
 552   erValue     = Interpolate( &fgkZList[fILow], savedEr, order, z )    ;
 553   ephiValue   = Interpolate( &fgkZList[fILow], savedEphi, order, z )  ;
 554   ezValue     = Interpolate( &fgkZList[fILow], savedEz, order, z )    ;
 555
 556 }
 557
 558 Double_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
 559                                               const Int_t nx,  const Int_t ny, const Double_t xv[], const Double_t yv[],
 560                                               const TMatrixD &array ) {
 561   //
 562   // Interpolate table (TMatrix format) - 2D interpolation
 563   //
 564
 565   static  Int_t jlow = 0, klow = 0 ;
 566   Double_t saveArray[10]  ;
 567
 568   Search( nx,  xv,  x,   jlow  ) ;
 569   Search( ny,  yv,  y,   klow  ) ;
 570   if ( jlow < 0 ) jlow = 0 ;   // check if out of range
 571   if ( klow < 0 ) klow = 0 ;
 572   if ( jlow + order  >=    nx - 1 ) jlow =   nx - 1 - order ;
 573   if ( klow + order  >=    ny - 1 ) klow =   ny - 1 - order ;
 574
 575   for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
 576     {
 577       Double_t *ajkl = &((TMatrixD&)array)(j,klow);
 578       saveArray[j-jlow]  = Interpolate( &yv[klow], ajkl , order, y )   ;
 579     }
 580
 581   return( Interpolate( &xv[jlow], saveArray, order, x ) )   ;
 582
 583 }
 584
 585 Double_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x,   const Double_t y,   const Double_t z,
 586                                               const Int_t  nx,    const Int_t  ny,    const Int_t  nz,
 587                                               const Double_t xv[], const Double_t yv[], const Double_t zv[],
 588                                               TMatrixD **arrayofArrays ) {
 589   //
 590   // Interpolate table (TMatrix format) - 3D interpolation
 591   //
 592
 593   static  Int_t ilow = 0, jlow = 0, klow = 0 ;
 594   Double_t saveArray[10],  savedArray[10] ;
 595
 596   Search( nx, xv, x, ilow   ) ;
 597   Search( ny, yv, y, jlow   ) ;
 598   Search( nz, zv, z, klow   ) ;
 599
 600   if ( ilow < 0 ) ilow = 0 ;   // check if out of range
 601   if ( jlow < 0 ) jlow = 0 ;
 602   if ( klow < 0 ) klow = 0 ;
 603
 604   if ( ilow + order  >=    nx - 1 ) ilow =   nx - 1 - order ;
 605   if ( jlow + order  >=    ny - 1 ) jlow =   ny - 1 - order ;
 606   if ( klow + order  >=    nz - 1 ) klow =   nz - 1 - order ;
 607
 608   for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
 609     {
 610       TMatrixD &table = *arrayofArrays[k] ;
 611       for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
 612         {
 613           saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y )   ;
 614         }
 615       savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x )   ;
 616     }
 617   return( Interpolate( &zv[klow], savedArray, order, z ) )   ;
 618
 619 }
 620
 621
 622
 623 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
 624                                        const Int_t order, const Double_t x ) {
 625   //
 626   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
 627   //
 628
 629   Double_t y ;
 630   if ( order == 2 ) {                // Quadratic Interpolation = 2
 631     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
 632     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
 633     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
 634   } else {                           // Linear Interpolation = 1
 635     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
 636   }
 637
 638   return (y);
 639
 640 }
 641
 642
 643 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
 644   //
 645   // Search an ordered table by starting at the most recently used point
 646   //
 647
 648   Long_t middle, high ;
 649   Int_t  ascend = 0, increment = 1 ;
 650
 651   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
 652
 653   if ( low < 0 || low > n-1 ) {
 654     low = -1 ; high = n ;
 655   } else {                                            // Ordered Search phase
 656     if ( (Int_t)( x >= xArray[low] ) == ascend )  {
 657       if ( low == n-1 ) return ;
 658       high = low + 1 ;
 659       while ( (Int_t)( x >= xArray[high] ) == ascend ) {
 660         low = high ;
 661         increment *= 2 ;
 662         high = low + increment ;
 663         if ( high > n-1 )  {  high = n ; break ;  }
 664       }
 665     } else {
 666       if ( low == 0 )  {  low = -1 ;  return ;  }
 667       high = low - 1 ;
 668       while ( (Int_t)( x < xArray[low] ) == ascend ) {
 669         high = low ;
 670         increment *= 2 ;
 671         if ( increment >= high )  {  low = -1 ;  break ;  }
 672         else  low = high - increment ;
 673       }
 674     }
 675   }
 676
 677   while ( (high-low) != 1 ) {                     // Binary Search Phase
 678     middle = ( high + low ) / 2 ;
 679     if ( (Int_t)( x >= xArray[middle] ) == ascend )
 680       low = middle ;
 681     else
 682       high = middle ;
 683   }
 684
 685   if ( x == xArray[n-1] ) low = n-2 ;
 686   if ( x == xArray[0]   ) low = 0 ;
 687
 688 }
 689
 690 void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
 691                                            TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
 692                                            const Int_t rows, const Int_t columns, const Int_t iterations,
 693                                            const Bool_t rocDisplacement ) {
 694   //
 695   // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
 696   //
 697   // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
 698   // boundary conditions on the first and last rows, and the first and last columns.  The remainder of the
 699   // array can be blank or contain a preliminary guess at the solution.  The Charge density matrix contains
 700   // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
 701   // you wish to solve Laplaces equation however it should not contain random numbers or you will get
 702   // random numbers back as a solution.
 703   // Poisson's equation is solved by iteratively relaxing the matrix to the final solution.  In order to
 704   // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
 705   // space.  First it solves the problem on a very sparse grid by skipping rows and columns in the original
 706   // matrix.  Then it doubles the number of points and solves the problem again.  Then it doubles the
 707   // number of points and solves the problem again.  This happens several times until the maximum number
 708   // of points has been included in the array.
 709   //
 710   // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
 711   // So rows == 2**M + 1 and columns == 2**N + 1.  The number of rows and columns can be different.
 712   //
 713   // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
 714   //
 715   // Original code by Jim Thomas (STAR TPC Collaboration)
 716   //
 717
 718   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
 719
 720   const Float_t  gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
 721   const Float_t  gridSizeZ   =  fgkTPCZ0 / (columns-1) ;
 722   const Float_t  ratio       =  gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
 723
 724   TMatrixD  arrayEr(rows,columns) ;
 725   TMatrixD  arrayEz(rows,columns) ;
 726
 727   //Check that number of rows and columns is suitable for a binary expansion
 728
 729   if ( !IsPowerOfTwo(rows-1) ) {
 730     AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
 731     return;
 732   }
 733   if ( !IsPowerOfTwo(columns-1) ) {
 734     AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
 735     return;
 736   }
 737
 738   // Solve Poisson's equation in cylindrical coordinates by relaxation technique
 739   // Allow for different size grid spacing in R and Z directions
 740   // Use a binary expansion of the size of the matrix to speed up the solution of the problem
 741
 742   Int_t iOne = (rows-1)/4 ;
 743   Int_t jOne = (columns-1)/4 ;
 744   // Solve for N in 2**N, add one.
 745   Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
 746
 747   for ( Int_t count = 0 ; count < loops ; count++ ) {
 748     // Loop while the matrix expands & the resolution increases.
 749
 750     Float_t tempGridSizeR = gridSizeR * iOne ;
 751     Float_t tempRatio     = ratio * iOne * iOne / ( jOne * jOne ) ;
 752     Float_t tempFourth    = 1.0 / (2.0 + 2.0*tempRatio) ;
 753
 754     // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
 755     std::vector<float> coef1(rows) ;
 756     std::vector<float> coef2(rows) ;
 757
 758     for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
 759        Float_t radius = fgkIFCRadius + i*gridSizeR ;
 760       coef1[i] = 1.0 + tempGridSizeR/(2*radius);
 761       coef2[i] = 1.0 - tempGridSizeR/(2*radius);
 762     }
 763
 764     TMatrixD sumChargeDensity(rows,columns) ;
 765
 766     for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 767       Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
 768       for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 769         if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
 770         else {
 771           // Add up all enclosed charge density contributions within 1/2 unit in all directions
 772           Float_t weight = 0.0 ;
 773           Float_t sum    = 0.0 ;
 774           sumChargeDensity(i,j) = 0.0 ;
 775           for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
 776             for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
 777               if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
 778               else
 779                 weight = 1.0 ;
 780               // Note that this is cylindrical geometry
 781               sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
 782               sum += weight*radius ;
 783             }
 784           }
 785           sumChargeDensity(i,j) /= sum ;
 786         }
 787         sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
 788        }
 789     }
 790
 791     for ( Int_t k = 1 ; k <= iterations; k++ ) {
 792       // Solve Poisson's Equation
 793       // Over-relaxation index, must be >= 1 but < 2.  Arrange for it to evolve from 2 => 1
 794       // as interations increase.
 795       Float_t overRelax   = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
 796       Float_t overRelaxM1 = overRelax - 1.0 ;
 797       Float_t overRelaxtempFourth, overRelaxcoef5 ;
 798       overRelaxtempFourth = overRelax * tempFourth ;
 799       overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
 800
 801       for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 802         for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 803
 804           arrayV(i,j) = (   coef2[i]       *   arrayV(i-iOne,j)
 805                           + tempRatio      * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
 806                           - overRelaxcoef5 *   arrayV(i,j)
 807                           + coef1[i]       *   arrayV(i+iOne,j)
 808                           + sumChargeDensity(i,j)
 809                         ) * overRelaxtempFourth;
 810         }
 811       }
 812
 813       if ( k == iterations ) {
 814         // After full solution is achieved, copy low resolution solution into higher res array
 815         for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 816           for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 817
 818             if ( iOne > 1 ) {
 819               arrayV(i+iOne/2,j)                    =  ( arrayV(i+iOne,j) + arrayV(i,j)     ) / 2 ;
 820               if ( i == iOne )  arrayV(i-iOne/2,j) =  ( arrayV(0,j)       + arrayV(iOne,j) ) / 2 ;
 821             }
 822             if ( jOne > 1 ) {
 823               arrayV(i,j+jOne/2)                    =  ( arrayV(i,j+jOne) + arrayV(i,j) )     / 2 ;
 824               if ( j == jOne )  arrayV(i,j-jOne/2) =  ( arrayV(i,0)       + arrayV(i,jOne) ) / 2 ;
 825             }
 826             if ( iOne > 1 && jOne > 1 ) {
 827               arrayV(i+iOne/2,j+jOne/2) =  ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
 828               if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
 829               if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
 830               // Note that this leaves a point at the upper left and lower right corners uninitialized.
 831               // -> Not a big deal.
 832             }
 833
 834           }
 835         }
 836       }
 837
 838     }
 839
 840     iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
 841     jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
 842
 843     sumChargeDensity.Clear();
 844   }
 845
 846   // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
 847   for ( Int_t j = 0 ; j < columns ; j++ ) {
 848     for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
 849     arrayEr(0,j)      =  -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
 850     arrayEr(rows-1,j) =  -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
 851   }
 852
 853   // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
 854   for ( Int_t i = 0 ; i < rows ; i++) {
 855     for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
 856     arrayEz(i,0)         =  -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
 857     arrayEz(i,columns-1) =  -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
 858   }
 859
 860   for ( Int_t i = 0 ; i < rows ; i++) {
 861     // Note: go back and compare to old version of this code.  See notes below.
 862     // JT Test ... attempt to divide by real Ez not Ez to first order
 863     for ( Int_t j = 0 ; j < columns ; j++ ) {
 864       arrayEz(i,j) += ezField;
 865       // This adds back the overall Z gradient of the field (main E field component)
 866     }
 867     // Warning: (-=) assumes you are using an error potetial without the overall Field included
 868   }
 869
 870   // Integrate Er/Ez from Z to zero
 871   for ( Int_t j = 0 ; j < columns ; j++ )  {
 872     for ( Int_t i = 0 ; i < rows ; i++ ) {
 873
 874       Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
 875       arrayErOverEz(i,j) = 0.0 ;
 876       arrayDeltaEz(i,j) = 0.0 ;
 877
 878       for ( Int_t k = j ; k < columns ; k++ ) {
 879         arrayErOverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
 880         arrayDeltaEz(i,j)   +=  index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
 881         if ( index != 4 )  index = 4; else index = 2 ;
 882       }
 883       if ( index == 4 ) {
 884         arrayErOverEz(i,j)  -=  (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
 885         arrayDeltaEz(i,j)   -=  (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
 886       }
 887       if ( index == 2 ) {
 888         arrayErOverEz(i,j)  +=  (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
 889                                                     -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
 890         arrayDeltaEz(i,j)   +=  (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
 891                                                     -2.5*(arrayEz(i,columns-1)-ezField));
 892       }
 893       if ( j == columns-2 ) {
 894         arrayErOverEz(i,j) =  (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
 895                                                   +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
 896         arrayDeltaEz(i,j)  =  (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
 897                                                   +1.5*(arrayEz(i,columns-1)-ezField) ) ;
 898       }
 899       if ( j == columns-1 ) {
 900         arrayErOverEz(i,j) =  0.0 ;
 901         arrayDeltaEz(i,j)  =  0.0 ;
 902       }
 903     }
 904   }
 905
 906   // calculate z distortion from the integrated Delta Ez residuals
 907   // and include the aquivalence (Volt to cm) of the ROC shift !!
 908
 909   for ( Int_t j = 0 ; j < columns ; j++ )  {
 910     for ( Int_t i = 0 ; i < rows ; i++ ) {
 911
 912       // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
 913       arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
 914
 915       // ROC Potential in cm aquivalent
 916       Double_t dzROCShift =  arrayV(i, columns -1)/ezField;
 917       if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift;  // add the ROC misaligment
 918
 919     }
 920   }
 921
 922   arrayEr.Clear();
 923   arrayEz.Clear();
 924
 925 }
 926
 927 void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
 928                     TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
 929                     const Int_t rows, const Int_t columns,  const Int_t phislices,
 930                     const Float_t deltaphi, const Int_t iterations, const Int_t symmetry,
 931                     Bool_t rocDisplacement  ) {
 932   //
 933   // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
 934   //
 935   //    NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
 936   //    The number of rows and COLUMNS can be different.
 937   //
 938   //    ROWS       ==  2**M + 1
 939   //    COLUMNS    ==  2**N + 1
 940   //    PHISLICES  ==  Arbitrary but greater than 3
 941   //
 942   //    DeltaPhi in Radians
 943   //
 944   //    SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
 945   //             = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
 946   //
 947   // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
 948
 949   const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
 950
 951   const Float_t  gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
 952   const Float_t  gridSizePhi =  deltaphi ;
 953   const Float_t  gridSizeZ   =  fgkTPCZ0 / (columns-1) ;
 954   const Float_t  ratioPhi    =  gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
 955   const Float_t  ratioZ      =  gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
 956
 957   TMatrixD arrayE(rows,columns) ;
 958
 959   // Check that the number of rows and columns is suitable for a binary expansion
 960   if ( !IsPowerOfTwo((rows-1))    ) {
 961     AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
 962     return; }
 963   if ( !IsPowerOfTwo((columns-1)) ) {
 964     AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
 965     return; }
 966   if ( phislices <= 3   )  {
 967     AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
 968     return; }
 969   if  ( phislices > 1000 ) {
 970     AliError("Poisson3D  phislices > 1000 is not allowed (nor wise) ");
 971     return; }
 972
 973   // Solve Poisson's equation in cylindrical coordinates by relaxation technique
 974   // Allow for different size grid spacing in R and Z directions
 975   // Use a binary expansion of the matrix to speed up the solution of the problem
 976
 977   Int_t loops, mplus, mminus, signplus, signminus  ;
 978   Int_t ione = (rows-1)/4 ;
 979   Int_t jone = (columns-1)/4 ;
 980   loops = TMath::Max(ione, jone) ;      // Calculate the number of loops for the binary expansion
 981   loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ;  // Solve for N in 2**N
 982
 983   TMatrixD* arrayofSumChargeDensities[1000] ;    // Create temporary arrays to store low resolution charge arrays
 984
 985   for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
 986
 987   for ( Int_t count = 0 ; count < loops ; count++ ) {      // START the master loop and do the binary expansion
 988
 989     Float_t  tempgridSizeR   =  gridSizeR  * ione ;
 990     Float_t  tempratioPhi    =  ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
 991     Float_t  tempratioZ      =  ratioZ   * ione * ione / ( jone * jone ) ;
 992
 993     std::vector<float> coef1(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
 994     std::vector<float> coef2(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
 995     std::vector<float> coef3(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
 996     std::vector<float> coef4(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
 997
 998     for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
 999       Float_t radius = fgkIFCRadius + i*gridSizeR ;
1000       coef1[i] = 1.0 + tempgridSizeR/(2*radius);
1001       coef2[i] = 1.0 - tempgridSizeR/(2*radius);
1002       coef3[i] = tempratioPhi/(radius*radius);
1003       coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
1004     }
1005
1006     for ( Int_t m = 0 ; m < phislices ; m++ ) {
1007       TMatrixD &chargeDensity    = *arrayofChargeDensities[m] ;
1008       TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
1009       for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
1010         Float_t radius = fgkIFCRadius + i*gridSizeR ;
1011         for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
1012           if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1013           else {           // Add up all enclosed charge density contributions within 1/2 unit in all directions
1014             Float_t weight = 0.0 ;
1015             Float_t sum    = 0.0 ;
1016             sumChargeDensity(i,j) = 0.0 ;
1017             for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
1018               for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
1019                 if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
1020                 else
1021                   weight = 1.0 ;
1022                 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1023                 sum += weight*radius ;
1024               }
1025             }
1026             sumChargeDensity(i,j) /= sum ;
1027           }
1028           sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
1029         }
1030       }
1031     }
1032
1033     for ( Int_t k = 1 ; k <= iterations; k++ ) {
1034
1035       // over-relaxation index, >= 1 but < 2
1036       Float_t overRelax   = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1037       Float_t overRelaxM1 = overRelax - 1.0 ;
1038
1039       std::vector<float> overRelaxcoef4(rows) ;  // Do this the standard C++ way to avoid gcc extensions
1040       std::vector<float> overRelaxcoef5(rows) ;  // Do this the standard C++ way to avoid gcc extensions
1041
1042       for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1043         overRelaxcoef4[i] = overRelax * coef4[i] ;
1044         overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
1045       }
1046
1047       for ( Int_t m = 0 ; m < phislices ; m++ ) {
1048
1049         mplus  = m + 1;   signplus  = 1 ;
1050         mminus = m - 1 ;  signminus = 1 ;
1051         if (symmetry==1) {  // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1052           if ( mplus  > phislices-1 ) mplus  = phislices - 2 ;
1053           if ( mminus < 0 )           mminus = 1 ;
1054         }
1055         else if (symmetry==-1) {   // Anti-symmetry in phi
1056           if ( mplus  > phislices-1 ) { mplus  = phislices - 2 ; signplus  = -1 ; }
1057           if ( mminus < 0 )           { mminus = 1 ;             signminus = -1 ; }
1058         }
1059                 else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1060           if ( mplus  > phislices-1 ) mplus  = m + 1 - phislices ;
1061           if ( mminus < 0 )           mminus = m - 1 + phislices ;
1062         }
1063         TMatrixD& arrayV    =  *arrayofArrayV[m] ;
1064         TMatrixD& arrayVP   =  *arrayofArrayV[mplus] ;
1065         TMatrixD& arrayVM   =  *arrayofArrayV[mminus] ;
1066         TMatrixD& sumChargeDensity =  *arrayofSumChargeDensities[m] ;
1067
1068         for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
1069           for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1070
1071             arrayV(i,j) = (   coef2[i]          *   arrayV(i-ione,j)
1072                             + tempratioZ        * ( arrayV(i,j-jone)  +  arrayV(i,j+jone) )
1073                             - overRelaxcoef5[i] *   arrayV(i,j)
1074                             + coef1[i]          *   arrayV(i+ione,j)
1075                             + coef3[i]          * ( signplus*arrayVP(i,j)       +  signminus*arrayVM(i,j) )
1076                             + sumChargeDensity(i,j)
1077                           ) * overRelaxcoef4[i] ;
1078             // Note: over-relax the solution at each step.  This speeds up the convergance.
1079
1080           }
1081         }
1082
1083         if ( k == iterations ) {   // After full solution is achieved, copy low resolution solution into higher res array
1084           for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
1085             for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1086
1087               if ( ione > 1 ) {
1088                 arrayV(i+ione/2,j)                    =  ( arrayV(i+ione,j) + arrayV(i,j)     ) / 2 ;
1089                 if ( i == ione )  arrayV(i-ione/2,j) =  ( arrayV(0,j)       + arrayV(ione,j) ) / 2 ;
1090               }
1091               if ( jone > 1 ) {
1092                 arrayV(i,j+jone/2)                    =  ( arrayV(i,j+jone) + arrayV(i,j) )     / 2 ;
1093                 if ( j == jone )  arrayV(i,j-jone/2) =  ( arrayV(i,0)       + arrayV(i,jone) ) / 2 ;
1094               }
1095               if ( ione > 1 && jone > 1 ) {
1096                 arrayV(i+ione/2,j+jone/2) =  ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
1097                 if ( i == ione ) arrayV(i-ione/2,j-jone/2) =   ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
1098                 if ( j == jone ) arrayV(i-ione/2,j-jone/2) =   ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
1099                 // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
1100               }
1101             }
1102           }
1103         }
1104
1105       }
1106     }
1107
1108     ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
1109     jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
1110
1111   }
1112
1113   //Differentiate V(r) and solve for E(r) using special equations for the first and last row
1114   //Integrate E(r)/E(z) from point of origin to pad plane
1115
1116   for ( Int_t m = 0 ; m < phislices ; m++ ) {
1117     TMatrixD& arrayV    =  *arrayofArrayV[m] ;
1118     TMatrixD& eroverEz  =  *arrayofEroverEz[m] ;
1119
1120     for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {  // Count backwards to facilitate integration over Z
1121
1122       // Differentiate in R
1123       for ( Int_t i = 1 ; i < rows-1 ; i++ )  arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1124       arrayE(0,j)      =  -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1125       arrayE(rows-1,j) =  -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1126       // Integrate over Z
1127       for ( Int_t i = 0 ; i < rows ; i++ ) {
1128         Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
1129         eroverEz(i,j) = 0.0 ;
1130         for ( Int_t k = j ; k < columns ; k++ ) {
1131
1132           eroverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1133           if ( index != 4 )  index = 4; else index = 2 ;
1134         }
1135         if ( index == 4 ) eroverEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1136         if ( index == 2 ) eroverEz(i,j)  +=
1137           (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1138         if ( j == columns-2 ) eroverEz(i,j) =
1139           (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1140         if ( j == columns-1 ) eroverEz(i,j) =  0.0 ;
1141       }
1142     }
1143     // if ( m == 0 ) { TCanvas*  c1 =  new TCanvas("erOverEz","erOverEz",50,50,840,600) ;  c1 -> cd() ;
1144     // eroverEz.Draw("surf") ; } // JT test
1145   }
1146
1147   //Differentiate V(r) and solve for E(phi)
1148   //Integrate E(phi)/E(z) from point of origin to pad plane
1149
1150   for ( Int_t m = 0 ; m < phislices ; m++ ) {
1151
1152     mplus  = m + 1;   signplus  = 1 ;
1153     mminus = m - 1 ;  signminus = 1 ;
1154     if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1155       if ( mplus  > phislices-1 ) mplus  = phislices - 2 ;
1156       if ( mminus < 0 )           mminus = 1 ;
1157     }
1158     else if (symmetry==-1) {       // Anti-symmetry in phi
1159       if ( mplus  > phislices-1 ) { mplus  = phislices - 2 ;  signplus  = -1 ; }
1160       if ( mminus < 0 )           { mminus = 1 ;                    signminus = -1 ; }
1161     }
1162     else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
1163       if ( mplus  > phislices-1 ) mplus  = m + 1 - phislices ;
1164       if ( mminus < 0 )           mminus = m - 1 + phislices ;
1165     }
1166     TMatrixD &arrayVP     =  *arrayofArrayV[mplus] ;
1167     TMatrixD &arrayVM     =  *arrayofArrayV[mminus] ;
1168     TMatrixD &ePhioverEz  =  *arrayofEPhioverEz[m] ;
1169     for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1170       // Differentiate in Phi
1171       for ( Int_t i = 0 ; i < rows ; i++ ) {
1172         Float_t radius = fgkIFCRadius + i*gridSizeR ;
1173         arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
1174       }
1175       // Integrate over Z
1176       for ( Int_t i = 0 ; i < rows ; i++ ) {
1177         Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
1178         ePhioverEz(i,j) = 0.0 ;
1179         for ( Int_t k = j ; k < columns ; k++ ) {
1180
1181           ePhioverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1182           if ( index != 4 )  index = 4; else index = 2 ;
1183         }
1184         if ( index == 4 ) ePhioverEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1185         if ( index == 2 ) ePhioverEz(i,j)  +=
1186           (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1187         if ( j == columns-2 ) ePhioverEz(i,j) =
1188           (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1189         if ( j == columns-1 ) ePhioverEz(i,j) =  0.0 ;
1190       }
1191     }
1192     // if ( m == 5 ) { TCanvas* c2 =  new TCanvas("arrayE","arrayE",50,50,840,600) ;  c2 -> cd() ;
1193     // arrayE.Draw("surf") ; } // JT test
1194   }
1195
1196
1197   // Differentiate V(r) and solve for E(z) using special equations for the first and last row
1198   // Integrate (E(z)-Ezstd) from point of origin to pad plane
1199
1200   for ( Int_t m = 0 ; m < phislices ; m++ ) {
1201     TMatrixD& arrayV   =  *arrayofArrayV[m] ;
1202     TMatrixD& deltaEz  =  *arrayofDeltaEz[m] ;
1203
1204     // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1205     for ( Int_t i = 0 ; i < rows ; i++) {
1206       for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1207       arrayE(i,0)         =  -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1208       arrayE(i,columns-1) =  -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1209     }
1210
1211     for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {  // Count backwards to facilitate integration over Z
1212       // Integrate over Z
1213       for ( Int_t i = 0 ; i < rows ; i++ ) {
1214         Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
1215         deltaEz(i,j) = 0.0 ;
1216         for ( Int_t k = j ; k < columns ; k++ ) {
1217           deltaEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k) ;
1218           if ( index != 4 )  index = 4; else index = 2 ;
1219         }
1220         if ( index == 4 ) deltaEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1) ;
1221         if ( index == 2 ) deltaEz(i,j)  +=
1222           (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
1223         if ( j == columns-2 ) deltaEz(i,j) =
1224           (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
1225         if ( j == columns-1 ) deltaEz(i,j) =  0.0 ;
1226       }
1227     }
1228     // if ( m == 0 ) { TCanvas*  c1 =  new TCanvas("erOverEz","erOverEz",50,50,840,600) ;  c1 -> cd() ;
1229     // eroverEz.Draw("surf") ; } // JT test
1230
1231     // calculate z distortion from the integrated Delta Ez residuals
1232     // and include the aquivalence (Volt to cm) of the ROC shift !!
1233
1234     for ( Int_t j = 0 ; j < columns ; j++ )  {
1235       for ( Int_t i = 0 ; i < rows ; i++ ) {
1236
1237         // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1238         deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
1239
1240         // ROC Potential in cm aquivalent
1241         Double_t dzROCShift =  arrayV(i, columns -1)/ezField;
1242         if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift;  // add the ROC misaligment
1243
1244       }
1245     }
1246
1247   } // end loop over phi
1248
1249
1250
1251   for ( Int_t k = 0 ; k < phislices ; k++ )
1252     {
1253       arrayofSumChargeDensities[k]->Delete() ;
1254     }
1255
1256
1257
1258   arrayE.Clear();
1259 }
1260
1261
1262 Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1263   //
1264   // Helperfunction: Check if integer is a power of 2
1265   //
1266   Int_t j = 0;
1267   while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
1268   if ( j == 1 ) return(1) ;  // True
1269   return(0) ;                // False
1270 }
1271
1272
1273 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
1274   //
1275   // Fit the track parameters - without and with distortion
1276   // 1. Space points in the TPC are simulated along the trajectory
1277   // 2. Space points distorted
1278   // 3. Fits  the non distorted and distroted track to the reference plane at refX
1279   // 4. For visualization and debugging  purposes the space points and tracks can be stored  in the tree - using the TTreeSRedirector functionality
1280   //
1281   // trackIn   - input track parameters
1282   // refX     - reference X to fit the track
1283   // dir      - direction - out=1 or in=-1
1284   // pcstream -  debug streamer to check the results
1285   //
1286   // see AliExternalTrackParam.h documentation:
1287   // track1.fP[0] - local y (rphi)
1288   // track1.fP[1] - z
1289   // track1.fP[2] - sinus of local inclination angle
1290   // track1.fP[3] - tangent of deep angle
1291   // track1.fP[4] - 1/pt
1292
1293   AliTPCROC * roc = AliTPCROC::Instance();
1294   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
1295   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
1296   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1297   const Double_t kMaxSnp = 0.85;
1298   const Double_t kSigmaY=0.1;
1299   const Double_t kSigmaZ=0.1;
1300   const Double_t kMaxR=500;
1301   const Double_t kMaxZ=500;
1302   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1303   Int_t npoints1=0;
1304   Int_t npoints2=0;
1305
1306   AliExternalTrackParam  track(trackIn); //
1307   // generate points
1308   AliTrackPointArray pointArray0(npoints0);
1309   AliTrackPointArray pointArray1(npoints0);
1310   Double_t xyz[3];
1311   if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
1312   //
1313   // simulate the track
1314   Int_t npoints=0;
1315   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
1316   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
1317     if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
1318     track.GetXYZ(xyz);
1319     xyz[0]+=gRandom->Gaus(0,0.00005);
1320     xyz[1]+=gRandom->Gaus(0,0.00005);
1321     xyz[2]+=gRandom->Gaus(0,0.00005);
1322     if (TMath::Abs(track.GetZ())>kMaxZ) break;
1323     if (TMath::Abs(track.GetX())>kMaxR) break;
1324     AliTrackPoint pIn0;                               // space point
1325     AliTrackPoint pIn1;
1326     Int_t sector= (xyz[2]>0)? 0:18;
1327     pointArray0.GetPoint(pIn0,npoints);
1328     pointArray1.GetPoint(pIn1,npoints);
1329     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
1330     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
1331     DistortPoint(distPoint, sector);
1332     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
1333     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
1334     //
1335     track.Rotate(alpha);
1336     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
1337     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
1338     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
1339     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
1340     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
1341     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
1342     pointArray0.AddPoint(npoints, &pIn0);
1343     pointArray1.AddPoint(npoints, &pIn1);
1344     npoints++;
1345     if (npoints>=npoints0) break;
1346   }
1347   if (npoints<npoints0/2) return 0;
1348   //
1349   // refit track
1350   //
1351   AliExternalTrackParam *track0=0;
1352   AliExternalTrackParam *track1=0;
1353   AliTrackPoint   point1,point2,point3;
1354   if (dir==1) {  //make seed inner
1355     pointArray0.GetPoint(point1,1);
1356     pointArray0.GetPoint(point2,30);
1357     pointArray0.GetPoint(point3,60);
1358   }
1359   if (dir==-1){ //make seed outer
1360     pointArray0.GetPoint(point1,npoints-60);
1361     pointArray0.GetPoint(point2,npoints-30);
1362     pointArray0.GetPoint(point3,npoints-1);
1363   }
1364   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
1365   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
1366
1367   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
1368     Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
1369     //
1370     AliTrackPoint pIn0;
1371     AliTrackPoint pIn1;
1372     pointArray0.GetPoint(pIn0,ipoint);
1373     pointArray1.GetPoint(pIn1,ipoint);
1374     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
1375     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
1376     //
1377     if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1378     if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1379     if (TMath::Abs(track0->GetZ())>kMaxZ) break;
1380     if (TMath::Abs(track0->GetX())>kMaxR) break;
1381     if (TMath::Abs(track1->GetZ())>kMaxZ) break;
1382     if (TMath::Abs(track1->GetX())>kMaxR) break;
1383
1384     track.GetXYZ(xyz);  // distorted track also propagated to the same reference radius
1385     //
1386     Double_t pointPos[2]={0,0};
1387     Double_t pointCov[3]={0,0,0};
1388     pointPos[0]=prot0.GetY();//local y
1389     pointPos[1]=prot0.GetZ();//local z
1390     pointCov[0]=prot0.GetCov()[3];//simay^2
1391     pointCov[1]=prot0.GetCov()[4];//sigmayz
1392     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
1393     if (!track0->Update(pointPos,pointCov)) break;
1394     //
1395     Double_t deltaX=prot1.GetX()-prot0.GetX();   // delta X
1396     Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp()));  // deltaY due  delta X
1397     Double_t deltaZX=deltaX*track1->GetTgl();                           // deltaZ due  delta X
1398
1399     pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
1400     pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
1401     pointCov[0]=prot1.GetCov()[3];//simay^2
1402     pointCov[1]=prot1.GetCov()[4];//sigmayz
1403     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
1404     if (!track1->Update(pointPos,pointCov)) break;
1405     npoints1++;
1406     npoints2++;
1407   }
1408   if (npoints2<npoints)  return 0;
1409   AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
1410   track1->Rotate(track0->GetAlpha());
1411   AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
1412
1413   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
1414     "point0.="<<&pointArray0<<   //  points
1415     "point1.="<<&pointArray1<<   //  distorted points
1416     "trackIn.="<<&track<<       //  original track
1417     "track0.="<<track0<<         //  fitted track
1418     "track1.="<<track1<<         //  fitted distorted track
1419     "\n";
1420   new(&trackIn) AliExternalTrackParam(*track0);
1421   delete track0;
1422   return track1;
1423 }
1424
1425
1426
1427
1428
1429 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
1430   //
1431   // create the distortion tree on a mesh with granularity given by step
1432   // return the tree with distortions at given position
1433   // Map is created on the mesh with given step size
1434   //
1435   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
1436   Float_t xyz[3];
1437   for (Double_t x= -250; x<250; x+=step){
1438     for (Double_t y= -250; y<250; y+=step){
1439       Double_t r    = TMath::Sqrt(x*x+y*y);
1440       if (r<80) continue;
1441       if (r>250) continue;
1442       for (Double_t z= -250; z<250; z+=step){
1443         Int_t roc=(z>0)?0:18;
1444         xyz[0]=x;
1445         xyz[1]=y;
1446         xyz[2]=z;
1447         Double_t phi  = TMath::ATan2(y,x);
1448         DistortPoint(xyz,roc);
1449         Double_t r1    = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1450         Double_t phi1  = TMath::ATan2(xyz[1],xyz[0]);
1451         if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
1452         if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
1453         Double_t dx = xyz[0]-x;
1454         Double_t dy = xyz[1]-y;
1455         Double_t dz = xyz[2]-z;
1456         Double_t dr=r1-r;
1457         Double_t drphi=(phi1-phi)*r;
1458         (*pcstream)<<"distortion"<<
1459           "x="<<x<<           // original position
1460           "y="<<y<<
1461           "z="<<z<<
1462           "r="<<r<<
1463           "phi="<<phi<<
1464           "x1="<<xyz[0]<<      // distorted position
1465           "y1="<<xyz[1]<<
1466           "z1="<<xyz[2]<<
1467           "r1="<<r1<<
1468           "phi1="<<phi1<<
1469           //
1470           "dx="<<dx<<          // delta position
1471           "dy="<<dy<<
1472           "dz="<<dz<<
1473           "dr="<<dr<<
1474           "drphi="<<drphi<<
1475           "\n";
1476       }
1477     }
1478   }
1479   delete pcstream;
1480   TFile f(Form("correction%s.root",GetName()));
1481   TTree * tree = (TTree*)f.Get("distortion");
1482   TTree * tree2= tree->CopyTree("1");
1483   tree2->SetName(Form("dist%s",GetName()));
1484   tree2->SetDirectory(0);
1485   delete tree;
1486   return tree2;
1487 }
1488
1489
1490
1491
1492 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
1493   //
1494   // Make a fit tree:
1495   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1496   // calculates partial distortions
1497   // Partial distortion is stored in the resulting tree
1498   // Output is storred in the file distortion_<dettype>_<partype>.root
1499   // Partial  distortion is stored with the name given by correction name
1500   //
1501   //
1502   // Parameters of function:
1503   // input     - input tree
1504   // dtype     - distortion type 0 - ITSTPC,  1 -TPCTRD, 2 - TPCvertex
1505   // ppype     - parameter type
1506   // corrArray - array with partial corrections
1507   // step      - skipe entries  - if 1 all entries processed - it is slow
1508   // debug     0 if debug on also space points dumped - it is slow
1509
1510   const Double_t kMaxSnp = 0.85;
1511   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1512   //  const Double_t kB2C=-0.299792458e-3;
1513   const Int_t kMinEntries=50;
1514   Double_t phi,theta, snp, mean,rms, entries;
1515   tinput->SetBranchAddress("theta",&theta);
1516   tinput->SetBranchAddress("phi", &phi);
1517   tinput->SetBranchAddress("snp",&snp);
1518   tinput->SetBranchAddress("mean",&mean);
1519   tinput->SetBranchAddress("rms",&rms);
1520   tinput->SetBranchAddress("entries",&entries);
1521   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
1522   //
1523   Int_t nentries=tinput->GetEntries();
1524   Int_t ncorr=corrArray->GetEntries();
1525   Double_t corrections[100]={0}; //
1526   Double_t tPar[5];
1527   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1528   Double_t refX=0;
1529   Int_t dir=0;
1530   if (dtype==0) {refX=85.; dir=-1;}
1531   if (dtype==1) {refX=275.; dir=1;}
1532   if (dtype==2) {refX=85.; dir=-1;}
1533   if (dtype==3) {refX=360.; dir=-1;}
1534   //
1535   for (Int_t ientry=0; ientry<nentries; ientry+=step){
1536     tinput->GetEntry(ientry);
1537     if (TMath::Abs(snp)>kMaxSnp) continue;
1538     tPar[0]=0;
1539     tPar[1]=theta*refX;
1540     tPar[2]=snp;
1541     tPar[3]=theta;
1542     tPar[4]=(gRandom->Rndm()-0.5)*0.02;  // should be calculated - non equal to 0
1543     Double_t bz=AliTrackerBase::GetBz();
1544     if (refX>10. && TMath::Abs(bz)>0.1 )  tPar[4]=snp/(refX*bz*kB2C*2);
1545     tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
1546     AliExternalTrackParam track(refX,phi,tPar,cov);
1547     Double_t xyz[3];
1548     track.GetXYZ(xyz);
1549     Int_t id=0;
1550     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1551     if (ptype==4 &&bz<0) mean*=-1;  // interpret as curvature
1552     (*pcstream)<<"fit"<<
1553       "bz="<<bz<<         // magnetic filed used
1554       "dtype="<<dtype<<   // detector match type
1555       "ptype="<<ptype<<   // parameter type
1556       "theta="<<theta<<   // theta
1557       "phi="<<phi<<       // phi
1558       "snp="<<snp<<       // snp
1559       "mean="<<mean<<     // mean dist value
1560       "rms="<<rms<<       // rms
1561       "gx="<<xyz[0]<<         // global position at reference
1562       "gy="<<xyz[1]<<         // global position at reference
1563       "gz="<<xyz[2]<<         // global position at reference
1564       "dRrec="<<dRrec<<      // delta Radius in reconstruction
1565       "id="<<id<<             // track id
1566       "entries="<<entries;// number of entries in bin
1567     //
1568     for (Int_t icorr=0; icorr<ncorr; icorr++) {
1569       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1570       corrections[icorr]=0;
1571       if (entries>kMinEntries){
1572         AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1573         AliExternalTrackParam *trackOut = 0;
1574         if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1575         if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1576         if (dtype==0) {refX=85.; dir=-1;}
1577         if (dtype==1) {refX=275.; dir=1;}
1578         if (dtype==2) {refX=0; dir=-1;}
1579         if (dtype==3) {refX=360.; dir=-1;}
1580         //
1581         if (trackOut){
1582           AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1583           trackOut->Rotate(trackIn.GetAlpha());
1584           trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1585           //
1586           corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1587           delete trackOut;
1588         }else{
1589           corrections[icorr]=0;
1590         }
1591         if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature
1592       }
1593       Double_t dRdummy=0;
1594       (*pcstream)<<"fit"<<
1595         Form("%s=",corr->GetName())<<corrections[icorr]<<   // dump correction value
1596         Form("dR%s=",corr->GetName())<<dRdummy;   // dump dummy correction value not needed for tracks
1597                                                   // for points it is neccessary
1598     }
1599     (*pcstream)<<"fit"<<"\n";
1600   }
1601   delete pcstream;
1602 }
1603
1604
1605
1606 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1607   //
1608   // Make a laser fit tree for global minimization
1609   //
1610   const Double_t cutErrY=0.1;
1611   const Double_t cutErrZ=0.1;
1612   const Double_t kEpsilon=0.00000001;
1613   TVectorD *vecdY=0;
1614   TVectorD *vecdZ=0;
1615   TVectorD *veceY=0;
1616   TVectorD *veceZ=0;
1617   AliTPCLaserTrack *ltr=0;
1618   AliTPCLaserTrack::LoadTracks();
1619   tree->SetBranchAddress("dY.",&vecdY);
1620   tree->SetBranchAddress("dZ.",&vecdZ);
1621   tree->SetBranchAddress("eY.",&veceY);
1622   tree->SetBranchAddress("eZ.",&veceZ);
1623   tree->SetBranchAddress("LTr.",&ltr);
1624   Int_t entries= tree->GetEntries();
1625   TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1626   Double_t bz=AliTrackerBase::GetBz();
1627   //
1628
1629   for (Int_t ientry=0; ientry<entries; ientry++){
1630     tree->GetEntry(ientry);
1631     if (!ltr->GetVecGX()){
1632       ltr->UpdatePoints();
1633     }
1634     TVectorD * delta= (itype==0)? vecdY:vecdZ;
1635     TVectorD * err= (itype==0)? veceY:veceZ;
1636
1637     for (Int_t irow=0; irow<159; irow++){
1638       Int_t nentries = 1000;
1639       if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1640       if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1641       Int_t dtype=4;
1642       Double_t phi   =(*ltr->GetVecPhi())[irow];
1643       Double_t theta =ltr->GetTgl();
1644       Double_t mean=delta->GetMatrixArray()[irow];
1645       Double_t gx=0,gy=0,gz=0;
1646       Double_t snp = (*ltr->GetVecP2())[irow];
1647       Double_t rms = 0.1+err->GetMatrixArray()[irow];
1648       gx = (*ltr->GetVecGX())[irow];
1649       gy = (*ltr->GetVecGY())[irow];
1650       gz = (*ltr->GetVecGZ())[irow];
1651       Int_t bundle= ltr->GetBundle();
1652       Double_t dRrec=0;
1653       //
1654       // get delta R used in reconstruction
1655       AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
1656       AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1657       const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1658       Double_t xyz0[3]={gx,gy,gz};
1659       Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1660       //
1661       // old ExB correction
1662       //
1663       if(recoParam&&recoParam->GetUseExBCorrection()) {
1664         Double_t xyz1[3]={gx,gy,gz};
1665         calib->GetExB()->Correct(xyz0,xyz1);
1666         Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1667         dRrec=oldR-newR;
1668       }
1669       if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1670         Float_t xyz1[3]={gx,gy,gz};
1671         Int_t sector=(gz>0)?0:18;
1672         correction->CorrectPoint(xyz1, sector);
1673         Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1674         dRrec=oldR-newR;
1675       }
1676
1677
1678       (*pcstream)<<"fit"<<
1679         "bz="<<bz<<         // magnetic filed used
1680         "dtype="<<dtype<<   // detector match type
1681         "ptype="<<itype<<   // parameter type
1682         "theta="<<theta<<   // theta
1683         "phi="<<phi<<       // phi
1684         "snp="<<snp<<       // snp
1685         "mean="<<mean<<     // mean dist value
1686         "rms="<<rms<<       // rms
1687         "gx="<<gx<<         // global position
1688         "gy="<<gy<<         // global position
1689         "gz="<<gz<<         // global position
1690         "dRrec="<<dRrec<<      // delta Radius in reconstruction
1691         "id="<<bundle<<     //bundle
1692         "entries="<<nentries;// number of entries in bin
1693       //
1694       //
1695       Double_t ky = TMath::Tan(TMath::ASin(snp));
1696       Int_t ncorr = corrArray->GetEntries();
1697       Double_t r0   = TMath::Sqrt(gx*gx+gy*gy);
1698       Double_t phi0 = TMath::ATan2(gy,gx);
1699       Double_t distortions[1000]={0};
1700       Double_t distortionsR[1000]={0};
1701       for (Int_t icorr=0; icorr<ncorr; icorr++) {
1702         AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1703         Float_t distPoint[3]={gx,gy,gz};
1704         Int_t sector= (gz>0)? 0:18;
1705         if (r0>80){
1706           corr->DistortPoint(distPoint, sector);
1707         }
1708         // Double_t value=distPoint[2]-gz;
1709         if (itype==0){
1710           Double_t r1   = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1711           Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1712           Double_t drphi= r0*(phi1-phi0);
1713           Double_t dr   = r1-r0;
1714           distortions[icorr]  = drphi-ky*dr;
1715           distortionsR[icorr] = dr;
1716         }
1717         (*pcstream)<<"fit"<<
1718           Form("%s=",corr->GetName())<<distortions[icorr]<<    // dump correction value
1719           Form("dR%s=",corr->GetName())<<distortionsR[icorr];   // dump correction R  value
1720       }
1721       (*pcstream)<<"fit"<<"\n";
1722     }
1723   }
1724   delete pcstream;
1725 }
1726
1727
1728
1729 void   AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
1730   //
1731   // make a distortion map out ou fthe residual histogram
1732   // Results are written to the debug streamer - pcstream
1733   // Parameters:
1734   //   his0       - input (4D) residual histogram
1735   //   pcstream   - file to write the tree
1736   //   run        - run number
1737   // marian.ivanov@cern.ch
1738   const Int_t kMinEntries=50;
1739   Int_t nbins1=his0->GetAxis(1)->GetNbins();
1740   Int_t first1=his0->GetAxis(1)->GetFirst();
1741   Int_t last1 =his0->GetAxis(1)->GetLast();
1742   //
1743   Double_t bz=AliTrackerBase::GetBz();
1744   Int_t idim[4]={0,1,2,3};
1745   for (Int_t ibin1=first1; ibin1<last1; ibin1++){   //axis 1 - theta
1746     //
1747     his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1748     Double_t       x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1749     THnSparse * his1 = his0->Projection(4,idim);  // projected histogram according range1
1750     Int_t nbins3     = his1->GetAxis(3)->GetNbins();
1751     Int_t first3     = his1->GetAxis(3)->GetFirst();
1752     Int_t last3      = his1->GetAxis(3)->GetLast();
1753     //
1754
1755     for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){   // axis 3 - local angle
1756       his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1757       Double_t      x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1758       if (ibin3<first3) {
1759         his1->GetAxis(3)->SetRangeUser(-1,1);
1760         x3=0;
1761       }
1762       THnSparse * his3= his1->Projection(4,idim);         //projected histogram according selection 3
1763       Int_t nbins2    = his3->GetAxis(2)->GetNbins();
1764       Int_t first2    = his3->GetAxis(2)->GetFirst();
1765       Int_t last2     = his3->GetAxis(2)->GetLast();
1766       //
1767       for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1768         his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1769         Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1770         TH1 * hisDelta = his3->Projection(0);
1771         //
1772         Double_t entries = hisDelta->GetEntries();
1773         Double_t mean=0, rms=0;
1774         if (entries>kMinEntries){
1775           mean    = hisDelta->GetMean();
1776           rms = hisDelta->GetRMS();
1777         }
1778         (*pcstream)<<hname<<
1779           "run="<<run<<
1780           "bz="<<bz<<
1781           "theta="<<x1<<
1782           "phi="<<x2<<
1783           "snp="<<x3<<
1784           "entries="<<entries<<
1785           "mean="<<mean<<
1786           "rms="<<rms<<
1787           "\n";
1788         delete hisDelta;
1789         printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1790       }
1791       delete his3;
1792     }
1793     delete his1;
1794   }
1795 }
1796
1797
1798
1799
1800
1801 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1802   //
1803   // Store object in the OCDB
1804   // By default the object is stored in the current directory
1805   // default comment consit of user name and the date
1806   //
1807   TString ocdbStorage="";
1808   ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1809   AliCDBMetaData *metaData= new AliCDBMetaData();
1810   metaData->SetObjectClassName("AliTPCCorrection");
1811   metaData->SetResponsible("Marian Ivanov");
1812   metaData->SetBeamPeriod(1);
1813   metaData->SetAliRootVersion("05-25-01"); //root version
1814   TString userName=gSystem->GetFromPipe("echo $USER");
1815   TString date=gSystem->GetFromPipe("date");
1816
1817   if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1818   if (comment) metaData->SetComment(comment);
1819   AliCDBId* id1=NULL;
1820   id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1821   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1822   gStorage->Put(this, (*id1), metaData);
1823 }
1824
1825
1826 void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
1827   //
1828   // Fast method to simulate the influence of the given distortion on the vertex reconstruction
1829   //
1830
1831   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
1832   if (!magF) AliError("Magneticd field - not initialized");
1833   Double_t bz = magF->SolenoidField(); //field in kGauss
1834   printf("bz: %lf\n",bz);
1835   AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
1836
1837   TObjArray   aTrk;              // Original Track array of Aside
1838   TObjArray   daTrk;             // Distorted Track array of A side
1839   UShort_t    *aId = new UShort_t[nTracks];      // A side Track ID
1840   TObjArray   cTrk;
1841   TObjArray   dcTrk;
1842   UShort_t    *cId = new UShort_t [nTracks];
1843   Int_t id=0;
1844   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1845   TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
1846   fpt.SetParameters(7.24,0.120);
1847   fpt.SetNpx(10000);
1848   for(Int_t nt=0; nt<nTracks; nt++){
1849     Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
1850     Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
1851     Double_t pt = fpt.GetRandom(); // momentum for f1
1852     //   printf("phi %lf  eta %lf pt %lf\n",phi,eta,pt);
1853     Short_t sign=1;
1854     if(gRandom->Rndm() < 0.5){
1855       sign =1;
1856     }else{
1857       sign=-1;
1858     }
1859
1860     Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1861     Double_t pxyz[3];
1862     pxyz[0]=pt*TMath::Cos(phi);
1863     pxyz[1]=pt*TMath::Sin(phi);
1864     pxyz[2]=pt*TMath::Tan(theta);
1865     Double_t cv[21]={0};
1866     AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
1867
1868     Double_t refX=1.;
1869     Int_t dir=-1;
1870     AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir,  NULL);
1871     if (!td) continue;
1872     if (pcstream) (*pcstream)<<"track"<<
1873       "eta="<<eta<<
1874       "theta="<<theta<<
1875       "tOrig.="<<t<<
1876       "td.="<<td<<
1877       "\n";
1878     if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
1879       if (td){
1880         daTrk.AddLast(td);
1881         aTrk.AddLast(t);
1882         Int_t nn=aTrk.GetEntriesFast();
1883         aId[nn]=id;
1884       }
1885     }else if(( eta<-0.07 )&&( eta>-etaCuts )){
1886       if (td){
1887         dcTrk.AddLast(td);
1888         cTrk.AddLast(t);
1889         Int_t nn=cTrk.GetEntriesFast();
1890         cId[nn]=id;
1891       }
1892     }
1893     id++;
1894   }// end of track loop
1895
1896   vertexer->SetTPCMode();
1897   vertexer->SetConstraintOff();
1898
1899   aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));
1900   avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
1901   cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));
1902   cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
1903   if (pcstream) (*pcstream)<<"vertex"<<
1904     "x="<<orgVertex[0]<<
1905     "y="<<orgVertex[1]<<
1906     "z="<<orgVertex[2]<<
1907     "av.="<<&aV<<              // distorted vertex A side
1908     "cv.="<<&cV<<              // distroted vertex C side
1909     "avO.="<<&avOrg<<         // original vertex A side
1910     "cvO.="<<&cvOrg<<
1911     "\n";
1912   delete []aId;
1913   delete []cId;
1914 }
1915
1916 void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
1917   //
1918   // make correction available for visualization using
1919   // TFormula, TFX and TTree::Draw
1920   // important in order to check corrections and also compute dervied variables
1921   // e.g correction partial derivatives
1922   //
1923   // NOTE - class is not owner of correction
1924   //
1925   if (!fgVisualCorrection) fgVisualCorrection=new TObjArray;
1926   if (position>=fgVisualCorrection->GetEntriesFast()) fgVisualCorrection->Expand(position*2);
1927   fgVisualCorrection->AddAt(corr, position);
1928 }
1929
1930
1931
1932 Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t r, Double_t kZ, Int_t axisType, Int_t corrType){
1933   //
1934   // calculate the correction at given position - check the geffCorr
1935   //
1936   if (!fgVisualCorrection) return 0;
1937   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1938   if (!corr) return 0;
1939   Double_t phi=sector*TMath::Pi()/9.;
1940   Double_t gx = r*TMath::Cos(phi);
1941   Double_t gy = r*TMath::Sin(phi);
1942   Double_t gz = r*kZ;
1943   Int_t nsector=(gz>0) ? 0:18;
1944   //
1945   //
1946   //
1947   Float_t distPoint[3]={gx,gy,gz};
1948   corr->DistortPoint(distPoint, nsector);
1949   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1950   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1951   Double_t phi0=TMath::ATan2(gy,gx);
1952   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1953   if (axisType==0) return r1-r0;
1954   if (axisType==1) return (phi1-phi0)*r0;
1955   if (axisType==2) return distPoint[2]-gz;
1956   return phi1-phi0;
1957 }
1958
1959 Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
1960   //
1961   // return correction at given x,y,z
1962   //
1963   if (!fgVisualCorrection) return 0;
1964   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1965   if (!corr) return 0;
1966   Double_t phi0= TMath::ATan2(gy,gx);
1967   Int_t nsector=(gz>0) ? 0:18;
1968   Float_t distPoint[3]={gx,gy,gz};
1969   corr->DistortPoint(distPoint, nsector);
1970   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1971   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1972   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1973   if (axisType==0) return r1-r0;
1974   if (axisType==1) return (phi1-phi0)*r0;
1975   if (axisType==2) return distPoint[2]-gz;
1976   return phi1-phi0;
1977 }