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