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
  55
  56 #include "TRandom.h"
  57 #include "AliExternalTrackParam.h"
  58 #include "AliTrackPointArray.h"
  59 #include "TDatabasePDG.h"
  60 #include "AliTrackerBase.h"
  61 #include "AliTPCROC.h"
  62 #include "THnSparse.h"
  63
  64 #include "TRandom.h"
  65 #include "AliTPCTransform.h"
  66 #include "AliTPCcalibDB.h"
  67 #include "AliTPCExB.h"
  68 #include "AliTPCCorrection.h"
  69 #include "AliTPCRecoParam.h"
  70
  71 #include  "AliExternalTrackParam.h"
  72 #include  "AliTrackPointArray.h"
  73 #include  "TDatabasePDG.h"
  74 #include  "AliTrackerBase.h"
  75 #include  "AliTPCROC.h"
  76 #include  "THnSparse.h"
  77
  78 #include  "AliTPCLaserTrack.h"
  79 #include  "AliESDVertex.h"
  80 #include  "AliVertexerTracks.h"
  81 #include  "TDatabasePDG.h"
  82 #include  "TF1.h"
  83
  84 #include "AliTPCCorrection.h"
  85 #include "AliLog.h"
  86
  87 ClassImp(AliTPCCorrection)
  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
  98 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
  99 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] =  {
 100 84.0,   84.5,   85.0,   85.5,   86.0,  87.0,    88.0,
 101 90.0,   92.0,   94.0,   96.0,   98.0,  100.0,  102.0,  104.0,  106.0,  108.0,
 102 110.0,  112.0,  114.0,  116.0,  118.0,  120.0,  122.0,  124.0,  126.0,  128.0,
 103 130.0,  132.0,  134.0,  136.0,  138.0,  140.0,  142.0,  144.0,  146.0,  148.0,
 104 150.0,  152.0,  154.0,  156.0,  158.0,  160.0,  162.0,  164.0,  166.0,  168.0,
 105 170.0,  172.0,  174.0,  176.0,  178.0,  180.0,  182.0,  184.0,  186.0,  188.0,
 106 190.0,  192.0,  194.0,  196.0,  198.0,  200.0,  202.0,  204.0,  206.0,  208.0,
 107 210.0,  212.0,  214.0,  216.0,  218.0,  220.0,  222.0,  224.0,  226.0,  228.0,
 108 230.0,  232.0,  234.0,  236.0,  238.0,  240.0,  242.0,  244.0,  246.0,  248.0,
 109 249.0,  249.5,  250.0,  251.5,  252.0  } ;
 110
 111 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ]     =   {
 112 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
 113 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
 114 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
 115 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
 116 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
 117 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
 118 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
 119 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
 120 -100.0,  -98.0,  -96.0,  -94.0,  -92.0,  -90.0,  -88.0,  -86.0,  -84.0,  -82.0,
 121 -80.0,  -78.0,  -76.0,  -74.0,  -72.0,  -70.0,  -68.0,  -66.0,  -64.0,  -62.0,
 122 -60.0,  -58.0,  -56.0,  -54.0,  -52.0,  -50.0,  -48.0,  -46.0,  -44.0,  -42.0,
 123 -40.0,  -38.0,  -36.0,  -34.0,  -32.0,  -30.0,  -28.0,  -26.0,  -24.0,  -22.0,
 124 -20.0,  -18.0,  -16.0,  -14.0,  -12.0,  -10.0,   -8.0,   -6.0,   -4.0,   -2.0,
 125 -1.0,   -0.5,   -0.2,   -0.1,  -0.05,   0.05,    0.1,    0.2,    0.5,    1.0,
 126  2.0,    4.0,    6.0,    8.0,   10.0,   12.0,   14.0,   16.0,   18.0,   20.0,
 127  22.0,   24.0,   26.0,   28.0,   30.0,   32.0,   34.0,   36.0,   38.0,   40.0,
 128  42.0,   44.0,   46.0,   48.0,   50.0,   52.0,   54.0,   56.0,   58.0,   60.0,
 129  62.0,   64.0,   66.0,   68.0,   70.0,   72.0,   74.0,   76.0,   78.0,   80.0,
 130  82.0,   84.0,   86.0,   88.0,   90.0,   92.0,   94.0,   96.0,   98.0,  100.0,
 131 102.0,  104.0,  106.0,  108.0,  110.0,  112.0,  114.0,  116.0,  118.0,  120.0,
 132 122.0,  124.0,  126.0,  128.0,  130.0,  132.0,  134.0,  136.0,  138.0,  140.0,
 133 142.0,  144.0,  146.0,  148.0,  150.0,  152.0,  154.0,  156.0,  158.0,  160.0,
 134 162.0,  164.0,  166.0,  168.0,  170.0,  172.0,  174.0,  176.0,  178.0,  180.0,
 135 182.0,  184.0,  186.0,  188.0,  190.0,  192.0,  194.0,  196.0,  198.0,  200.0,
 136 202.0,  204.0,  206.0,  208.0,  210.0,  212.0,  214.0,  216.0,  218.0,  220.0,
 137 222.0,  224.0,  226.0,  228.0,  230.0,  232.0,  234.0,  236.0,  238.0,  240.0,
 138 242.0,  243.0,  244.0,  245.0,  246.0,  247.0,  248.0,  248.5,  249.0,  249.5   } ;
 139
 140
 141
 142 AliTPCCorrection::AliTPCCorrection()
 143   : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
 144 {
 145   //
 146   // default constructor
 147   //
 148 }
 149
 150 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
 151 : TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
 152 {
 153   //
 154   // default constructor, that set the name and title
 155   //
 156 }
 157
 158 AliTPCCorrection::~AliTPCCorrection() {
 159   //
 160   // virtual destructor
 161   //
 162 }
 163
 164 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
 165   //
 166   // Corrects the initial coordinates x (cartesian coordinates)
 167   // according to the given effect (inherited classes)
 168   // roc represents the TPC read out chamber (offline numbering convention)
 169   //
 170   Float_t dx[3];
 171   GetCorrection(x,roc,dx);
 172   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 173 }
 174
 175 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 176   //
 177   // Corrects the initial coordinates x (cartesian coordinates) and stores the new
 178   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 179   // roc represents the TPC read out chamber (offline numbering convention)
 180   //
 181   Float_t dx[3];
 182   GetCorrection(x,roc,dx);
 183   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 184 }
 185
 186 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
 187   //
 188   // Distorts 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   GetDistortion(x,roc,dx);
 194   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 195 }
 196
 197 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 198   //
 199   // Distorts 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   GetDistortion(x,roc,dx);
 205   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 206 }
 207
 208 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
 209   //
 210   // This function delivers the correction values dx in respect to the inital coordinates x
 211   // roc represents the TPC read out chamber (offline numbering convention)
 212   // Note: The dx is overwritten by the inherited effectice class ...
 213   //
 214   for (Int_t j=0;j<3;++j) { dx[j]=0.; }
 215 }
 216
 217 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
 218   //
 219   // This function delivers the distortion values dx in respect to the inital coordinates x
 220   // roc represents the TPC read out chamber (offline numbering convention)
 221   //
 222   GetCorrection(x,roc,dx);
 223   for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
 224 }
 225
 226 void AliTPCCorrection::Init() {
 227   //
 228   // Initialization funtion (not used at the moment)
 229   //
 230 }
 231
 232 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
 233   //
 234   // Update function
 235   //
 236 }
 237
 238 void AliTPCCorrection::Print(Option_t* /*option*/) const {
 239   //
 240   // Print function to check which correction classes are used
 241   // option=="d" prints details regarding the setted magnitude
 242   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 243   //
 244   printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
 245 }
 246
 247 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
 248   //
 249   // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
 250   // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
 251   // calibration run
 252   //
 253   fT1=t1;
 254   fT2=t2;
 255   //SetOmegaTauT1T2(omegaTau, t1, t2);
 256 }
 257
 258 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
 259   //
 260   // Simple plot functionality.
 261   // Returns a 2d hisogram which represents the corrections in radial direction (dr)
 262   // in respect to position z within the XY plane.
 263   // The histogramm has nx times ny entries.
 264   //
 265
 266   TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
 267                      nx,-250.,250.,ny,-250.,250.);
 268   Float_t x[3],dx[3];
 269   x[2]=z;
 270   Int_t roc=z>0.?0:18; // FIXME
 271   for (Int_t iy=1;iy<=ny;++iy) {
 272     x[1]=h->GetYaxis()->GetBinCenter(iy);
 273     for (Int_t ix=1;ix<=nx;++ix) {
 274       x[0]=h->GetXaxis()->GetBinCenter(ix);
 275       GetCorrection(x,roc,dx);
 276       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 277       if (90.<=r0 && r0<=250.) {
 278         Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 279         h->SetBinContent(ix,iy,r1-r0);
 280       }
 281       else
 282         h->SetBinContent(ix,iy,0.);
 283     }
 284   }
 285   return h;
 286 }
 287
 288 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
 289   //
 290   // Simple plot functionality.
 291   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 292   // in respect to position z within the XY plane.
 293   // The histogramm has nx times ny entries.
 294   //
 295
 296   TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
 297                      nx,-250.,250.,ny,-250.,250.);
 298   Float_t x[3],dx[3];
 299   x[2]=z;
 300   Int_t roc=z>0.?0:18; // FIXME
 301   for (Int_t iy=1;iy<=ny;++iy) {
 302     x[1]=h->GetYaxis()->GetBinCenter(iy);
 303     for (Int_t ix=1;ix<=nx;++ix) {
 304       x[0]=h->GetXaxis()->GetBinCenter(ix);
 305       GetCorrection(x,roc,dx);
 306       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 307       if (90.<=r0 && r0<=250.) {
 308         Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 309         Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 310
 311         Float_t dphi=phi1-phi0;
 312         if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 313         if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 314
 315         h->SetBinContent(ix,iy,r0*dphi);
 316       }
 317       else
 318         h->SetBinContent(ix,iy,0.);
 319     }
 320   }
 321   return h;
 322 }
 323
 324 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
 325   //
 326   // Simple plot functionality.
 327   // Returns a 2d hisogram which represents the corrections in r direction (dr)
 328   // in respect to angle phi within the ZR plane.
 329   // The histogramm has nx times ny entries.
 330   //
 331   TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
 332                      nz,-250.,250.,nr,85.,250.);
 333   Float_t x[3],dx[3];
 334   for (Int_t ir=1;ir<=nr;++ir) {
 335     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 336     x[0]=radius*TMath::Cos(phi);
 337     x[1]=radius*TMath::Sin(phi);
 338     for (Int_t iz=1;iz<=nz;++iz) {
 339       x[2]=h->GetXaxis()->GetBinCenter(iz);
 340       Int_t roc=x[2]>0.?0:18; // FIXME
 341       GetCorrection(x,roc,dx);
 342       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 343       Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 344       h->SetBinContent(iz,ir,r1-r0);
 345     }
 346   }
 347   return h;
 348
 349 }
 350
 351 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
 352   //
 353   // Simple plot functionality.
 354   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 355   // in respect to angle phi within the ZR plane.
 356   // The histogramm has nx times ny entries.
 357   //
 358   TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
 359                      nz,-250.,250.,nr,85.,250.);
 360   Float_t x[3],dx[3];
 361   for (Int_t iz=1;iz<=nz;++iz) {
 362     x[2]=h->GetXaxis()->GetBinCenter(iz);
 363     Int_t roc=x[2]>0.?0:18; // FIXME
 364     for (Int_t ir=1;ir<=nr;++ir) {
 365       Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 366       x[0]=radius*TMath::Cos(phi);
 367       x[1]=radius*TMath::Sin(phi);
 368       GetCorrection(x,roc,dx);
 369       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 370       Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 371       Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 372
 373       Float_t dphi=phi1-phi0;
 374       if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 375       if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 376
 377       h->SetBinContent(iz,ir,r0*dphi);
 378     }
 379   }
 380   return h;
 381 }
 382
 383 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
 384                                    const char *xlabel,const char *ylabel,const char *zlabel,
 385                                   Int_t nbinsx,Double_t xlow,Double_t xup,
 386                                   Int_t nbinsy,Double_t ylow,Double_t yup) {
 387   //
 388   // Helper function to create a 2d histogramm of given size
 389   //
 390
 391   TString hname=name;
 392   Int_t i=0;
 393   if (gDirectory) {
 394     while (gDirectory->FindObject(hname.Data())) {
 395       hname =name;
 396       hname+="_";
 397       hname+=i;
 398       ++i;
 399     }
 400   }
 401   TH2F *h=new TH2F(hname.Data(),title,
 402                    nbinsx,xlow,xup,
 403                    nbinsy,ylow,yup);
 404   h->GetXaxis()->SetTitle(xlabel);
 405   h->GetYaxis()->SetTitle(ylabel);
 406   h->GetZaxis()->SetTitle(zlabel);
 407   h->SetStats(0);
 408   return h;
 409 }
 410
 411
 412 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
 413
 414 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
 415                                                   const Double_t er[kNZ][kNR], Double_t &erValue ) {
 416   //
 417   // Interpolate table - 2D interpolation
 418   //
 419   Double_t saveEr[10] ;
 420
 421   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
 422   Search( kNR,   fgkRList,  r,   fKLow   ) ;
 423   if ( fJLow < 0 ) fJLow = 0 ;   // check if out of range
 424   if ( fKLow < 0 ) fKLow = 0 ;
 425   if ( fJLow + order  >=    kNZ - 1 ) fJLow =   kNZ - 1 - order ;
 426   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 427
 428   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
 429       saveEr[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
 430   }
 431   erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z )   ;
 432
 433 }
 434
 435
 436 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
 437                                        const Int_t order, const Double_t x ) {
 438   //
 439   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
 440   //
 441
 442   Double_t y ;
 443   if ( order == 2 ) {                // Quadratic Interpolation = 2
 444     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
 445     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
 446     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
 447   } else {                           // Linear Interpolation = 1
 448     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
 449   }
 450
 451   return (y);
 452
 453 }
 454
 455
 456 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
 457   //
 458   // Search an ordered table by starting at the most recently used point
 459   //
 460
 461   Long_t middle, high ;
 462   Int_t  ascend = 0, increment = 1 ;
 463
 464   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
 465
 466   if ( low < 0 || low > n-1 ) {
 467     low = -1 ; high = n ;
 468   } else {                                            // Ordered Search phase
 469     if ( (Int_t)( x >= xArray[low] ) == ascend )  {
 470       if ( low == n-1 ) return ;
 471       high = low + 1 ;
 472       while ( (Int_t)( x >= xArray[high] ) == ascend ) {
 473         low = high ;
 474         increment *= 2 ;
 475         high = low + increment ;
 476         if ( high > n-1 )  {  high = n ; break ;  }
 477       }
 478     } else {
 479       if ( low == 0 )  {  low = -1 ;  return ;  }
 480       high = low - 1 ;
 481       while ( (Int_t)( x < xArray[low] ) == ascend ) {
 482         high = low ;
 483         increment *= 2 ;
 484         if ( increment >= high )  {  low = -1 ;  break ;  }
 485         else  low = high - increment ;
 486       }
 487     }
 488   }
 489
 490   while ( (high-low) != 1 ) {                     // Binary Search Phase
 491     middle = ( high + low ) / 2 ;
 492     if ( (Int_t)( x >= xArray[middle] ) == ascend )
 493       low = middle ;
 494     else
 495       high = middle ;
 496   }
 497
 498   if ( x == xArray[n-1] ) low = n-2 ;
 499   if ( x == xArray[0]   ) low = 0 ;
 500
 501 }
 502
 503 void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, const TMatrixD &chargeDensity,
 504                                            TMatrixD &arrayErOverEz, const Int_t rows,
 505                                            const Int_t columns, const Int_t iterations ) {
 506   //
 507   // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
 508   //
 509   // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
 510   // boundary conditions on the first and last rows, and the first and last columns.  The remainder of the
 511   // array can be blank or contain a preliminary guess at the solution.  The Charge density matrix contains
 512   // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
 513   // you wish to solve Laplaces equation however it should not contain random numbers or you will get
 514   // random numbers back as a solution.
 515   // Poisson's equation is solved by iteratively relaxing the matrix to the final solution.  In order to
 516   // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
 517   // space.  First it solves the problem on a very sparse grid by skipping rows and columns in the original
 518   // matrix.  Then it doubles the number of points and solves the problem again.  Then it doubles the
 519   // number of points and solves the problem again.  This happens several times until the maximum number
 520   // of points has been included in the array.
 521   //
 522   // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
 523   // So rows == 2**M + 1 and columns == 2**N + 1.  The number of rows and columns can be different.
 524   //
 525   // Original code by Jim Thomas (STAR TPC Collaboration)
 526   //
 527
 528   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
 529
 530   const Float_t  gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
 531   const Float_t  gridSizeZ   =  fgkTPCZ0 / (columns-1) ;
 532   const Float_t  ratio       =  gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
 533
 534   TMatrixD  arrayEr(rows,columns) ;
 535   TMatrixD  arrayEz(rows,columns) ;
 536
 537   //Check that number of rows and columns is suitable for a binary expansion
 538
 539   if ( !IsPowerOfTwo(rows-1) ) {
 540     AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
 541     return;
 542   }
 543   if ( !IsPowerOfTwo(columns-1) ) {
 544     AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
 545     return;
 546   }
 547
 548   // Solve Poisson's equation in cylindrical coordinates by relaxation technique
 549   // Allow for different size grid spacing in R and Z directions
 550   // Use a binary expansion of the size of the matrix to speed up the solution of the problem
 551
 552   Int_t iOne = (rows-1)/4 ;
 553   Int_t jOne = (columns-1)/4 ;
 554   // Solve for N in 2**N, add one.
 555   Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
 556
 557   for ( Int_t count = 0 ; count < loops ; count++ ) {
 558     // Loop while the matrix expands & the resolution increases.
 559
 560     Float_t tempGridSizeR = gridSizeR * iOne ;
 561     Float_t tempRatio     = ratio * iOne * iOne / ( jOne * jOne ) ;
 562     Float_t tempFourth    = 1.0 / (2.0 + 2.0*tempRatio) ;
 563
 564     // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
 565     std::vector<float> coef1(rows) ;
 566     std::vector<float> coef2(rows) ;
 567
 568     for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
 569        Float_t radius = fgkIFCRadius + i*gridSizeR ;
 570       coef1[i] = 1.0 + tempGridSizeR/(2*radius);
 571       coef2[i] = 1.0 - tempGridSizeR/(2*radius);
 572     }
 573
 574     TMatrixD sumChargeDensity(rows,columns) ;
 575
 576     for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 577       Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
 578       for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 579         if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
 580         else {
 581           // Add up all enclosed charge density contributions within 1/2 unit in all directions
 582           Float_t weight = 0.0 ;
 583           Float_t sum    = 0.0 ;
 584           sumChargeDensity(i,j) = 0.0 ;
 585           for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
 586             for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
 587               if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
 588               else
 589                 weight = 1.0 ;
 590               // Note that this is cylindrical geometry
 591               sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
 592               sum += weight*radius ;
 593             }
 594           }
 595           sumChargeDensity(i,j) /= sum ;
 596         }
 597         sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
 598        }
 599     }
 600
 601     for ( Int_t k = 1 ; k <= iterations; k++ ) {
 602       // Solve Poisson's Equation
 603       // Over-relaxation index, must be >= 1 but < 2.  Arrange for it to evolve from 2 => 1
 604       // as interations increase.
 605       Float_t overRelax   = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
 606       Float_t overRelaxM1 = overRelax - 1.0 ;
 607       Float_t overRelaxtempFourth, overRelaxcoef5 ;
 608       overRelaxtempFourth = overRelax * tempFourth ;
 609       overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
 610
 611       for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 612         for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 613
 614           arrayV(i,j) = (   coef2[i]       *   arrayV(i-iOne,j)
 615                           + tempRatio      * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
 616                           - overRelaxcoef5 *   arrayV(i,j)
 617                           + coef1[i]       *   arrayV(i+iOne,j)
 618                           + sumChargeDensity(i,j)
 619                         ) * overRelaxtempFourth;
 620         }
 621       }
 622
 623       if ( k == iterations ) {
 624         // After full solution is achieved, copy low resolution solution into higher res array
 625         for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
 626           for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
 627
 628             if ( iOne > 1 ) {
 629               arrayV(i+iOne/2,j)                    =  ( arrayV(i+iOne,j) + arrayV(i,j)     ) / 2 ;
 630               if ( i == iOne )  arrayV(i-iOne/2,j) =  ( arrayV(0,j)       + arrayV(iOne,j) ) / 2 ;
 631             }
 632             if ( jOne > 1 ) {
 633               arrayV(i,j+jOne/2)                    =  ( arrayV(i,j+jOne) + arrayV(i,j) )     / 2 ;
 634               if ( j == jOne )  arrayV(i,j-jOne/2) =  ( arrayV(i,0)       + arrayV(i,jOne) ) / 2 ;
 635             }
 636             if ( iOne > 1 && jOne > 1 ) {
 637               arrayV(i+iOne/2,j+jOne/2) =  ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
 638               if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
 639               if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
 640               // Note that this leaves a point at the upper left and lower right corners uninitialized.
 641               // -> Not a big deal.
 642             }
 643
 644           }
 645         }
 646       }
 647
 648     }
 649
 650     iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
 651     jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
 652
 653   }
 654
 655   // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
 656   for ( Int_t j = 0 ; j < columns ; j++ ) {
 657     for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
 658     arrayEr(0,j)      =  -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
 659     arrayEr(rows-1,j) =  -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
 660   }
 661
 662   // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
 663   for ( Int_t i = 0 ; i < rows ; i++) {
 664     for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
 665     arrayEz(i,0)         =  -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
 666     arrayEz(i,columns-1) =  -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
 667   }
 668
 669   for ( Int_t i = 0 ; i < rows ; i++) {
 670     // Note: go back and compare to old version of this code.  See notes below.
 671     // JT Test ... attempt to divide by real Ez not Ez to first order
 672     for ( Int_t j = 0 ; j < columns ; j++ ) {
 673       arrayEz(i,j) += ezField;
 674       // This adds back the overall Z gradient of the field (main E field component)
 675     }
 676     // Warning: (-=) assumes you are using an error potetial without the overall Field included
 677   }
 678
 679   // Integrate Er/Ez from Z to zero
 680   for ( Int_t j = 0 ; j < columns ; j++ )  {
 681     for ( Int_t i = 0 ; i < rows ; i++ ) {
 682       Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
 683       arrayErOverEz(i,j) = 0.0 ;
 684       for ( Int_t k = j ; k < columns ; k++ ) {
 685         arrayErOverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
 686         if ( index != 4 )  index = 4; else index = 2 ;
 687       }
 688       if ( index == 4 ) arrayErOverEz(i,j)  -=  (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
 689       if ( index == 2 ) arrayErOverEz(i,j)  +=
 690         (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
 691                             -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) )   ;
 692       if ( j == columns-2 ) arrayErOverEz(i,j) =
 693         (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
 694                             +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
 695       if ( j == columns-1 ) arrayErOverEz(i,j) =  0.0 ;
 696     }
 697   }
 698
 699 }
 700
 701
 702
 703 Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
 704   //
 705   // Helperfunction: Check if integer is a power of 2
 706   //
 707   Int_t j = 0;
 708   while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
 709   if ( j == 1 ) return(1) ;  // True
 710   return(0) ;                // False
 711 }
 712
 713
 714 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
 715   //
 716   // Fit the track parameters - without and with distortion
 717   // 1. Space points in the TPC are simulated along the trajectory
 718   // 2. Space points distorted
 719   // 3. Fits  the non distorted and distroted track to the reference plane at refX
 720   // 4. For visualization and debugging  purposes the space points and tracks can be stored  in the tree - using the TTreeSRedirector functionality
 721   //
 722   // trackIn   - input track parameters
 723   // refX     - reference X to fit the track
 724   // dir      - direction - out=1 or in=-1
 725   // pcstream -  debug streamer to check the results
 726   //
 727   // see AliExternalTrackParam.h documentation:
 728   // track1.fP[0] - local y (rphi)
 729   // track1.fP[1] - z
 730   // track1.fP[2] - sinus of local inclination angle
 731   // track1.fP[3] - tangent of deep angle
 732   // track1.fP[4] - 1/pt
 733
 734   AliTPCROC * roc = AliTPCROC::Instance();
 735   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
 736   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
 737   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 738   const Double_t kMaxSnp = 0.85;
 739   const Double_t kSigmaY=0.1;
 740   const Double_t kSigmaZ=0.1;
 741   const Double_t kMaxR=500;
 742   const Double_t kMaxZ=500;
 743   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 744   Int_t npoints1=0;
 745   Int_t npoints2=0;
 746
 747   AliExternalTrackParam  track(trackIn); //
 748   // generate points
 749   AliTrackPointArray pointArray0(npoints0);
 750   AliTrackPointArray pointArray1(npoints0);
 751   Double_t xyz[3];
 752   if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
 753   //
 754   // simulate the track
 755   Int_t npoints=0;
 756   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
 757   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
 758     if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
 759     track.GetXYZ(xyz);
 760     xyz[0]+=gRandom->Gaus(0,0.00005);
 761     xyz[1]+=gRandom->Gaus(0,0.00005);
 762     xyz[2]+=gRandom->Gaus(0,0.00005);
 763     if (TMath::Abs(track.GetZ())>kMaxZ) break;
 764     if (TMath::Abs(track.GetX())>kMaxR) break;
 765     AliTrackPoint pIn0;                               // space point
 766     AliTrackPoint pIn1;
 767     Int_t sector= (xyz[2]>0)? 0:18;
 768     pointArray0.GetPoint(pIn0,npoints);
 769     pointArray1.GetPoint(pIn1,npoints);
 770     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 771     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
 772     DistortPoint(distPoint, sector);
 773     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
 774     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
 775     //
 776     track.Rotate(alpha);
 777     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
 778     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
 779     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
 780     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
 781     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
 782     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
 783     pointArray0.AddPoint(npoints, &pIn0);
 784     pointArray1.AddPoint(npoints, &pIn1);
 785     npoints++;
 786     if (npoints>=npoints0) break;
 787   }
 788   if (npoints<npoints0/2) return 0;
 789   //
 790   // refit track
 791   //
 792   AliExternalTrackParam *track0=0;
 793   AliExternalTrackParam *track1=0;
 794   AliTrackPoint   point1,point2,point3;
 795   if (dir==1) {  //make seed inner
 796     pointArray0.GetPoint(point1,1);
 797     pointArray0.GetPoint(point2,30);
 798     pointArray0.GetPoint(point3,60);
 799   }
 800   if (dir==-1){ //make seed outer
 801     pointArray0.GetPoint(point1,npoints-60);
 802     pointArray0.GetPoint(point2,npoints-30);
 803     pointArray0.GetPoint(point3,npoints-1);
 804   }
 805   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
 806   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
 807
 808   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
 809     Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
 810     //
 811     AliTrackPoint pIn0;
 812     AliTrackPoint pIn1;
 813     pointArray0.GetPoint(pIn0,ipoint);
 814     pointArray1.GetPoint(pIn1,ipoint);
 815     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
 816     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
 817     //
 818     if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
 819     if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
 820     if (TMath::Abs(track0->GetZ())>kMaxZ) break;
 821     if (TMath::Abs(track0->GetX())>kMaxR) break;
 822     if (TMath::Abs(track1->GetZ())>kMaxZ) break;
 823     if (TMath::Abs(track1->GetX())>kMaxR) break;
 824
 825     track.GetXYZ(xyz);  // distorted track also propagated to the same reference radius
 826     //
 827     Double_t pointPos[2]={0,0};
 828     Double_t pointCov[3]={0,0,0};
 829     pointPos[0]=prot0.GetY();//local y
 830     pointPos[1]=prot0.GetZ();//local z
 831     pointCov[0]=prot0.GetCov()[3];//simay^2
 832     pointCov[1]=prot0.GetCov()[4];//sigmayz
 833     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
 834     if (!track0->Update(pointPos,pointCov)) break;
 835     //
 836     Double_t deltaX=prot1.GetX()-prot0.GetX();   // delta X
 837     Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp()));  // deltaY due  delta X
 838     Double_t deltaZX=deltaX*track1->GetTgl();                           // deltaZ due  delta X
 839
 840     pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
 841     pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
 842     pointCov[0]=prot1.GetCov()[3];//simay^2
 843     pointCov[1]=prot1.GetCov()[4];//sigmayz
 844     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
 845     if (!track1->Update(pointPos,pointCov)) break;
 846     npoints1++;
 847     npoints2++;
 848   }
 849   if (npoints2<npoints)  return 0;
 850   AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
 851   track1->Rotate(track0->GetAlpha());
 852   AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
 853
 854   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
 855     "point0.="<<&pointArray0<<   //  points
 856     "point1.="<<&pointArray1<<   //  distorted points
 857     "trackIn.="<<&track<<       //  original track
 858     "track0.="<<track0<<         //  fitted track
 859     "track1.="<<track1<<         //  fitted distorted track
 860     "\n";
 861   new(&trackIn) AliExternalTrackParam(*track0);
 862   delete track0;
 863   return track1;
 864 }
 865
 866
 867
 868
 869
 870 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
 871   //
 872   // create the distortion tree on a mesh with granularity given by step
 873   // return the tree with distortions at given position
 874   // Map is created on the mesh with given step size
 875   //
 876   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
 877   Float_t xyz[3];
 878   for (Double_t x= -250; x<250; x+=step){
 879     for (Double_t y= -250; y<250; y+=step){
 880       Double_t r    = TMath::Sqrt(x*x+y*y);
 881       if (r<80) continue;
 882       if (r>250) continue;
 883       for (Double_t z= -250; z<250; z+=step){
 884         Int_t roc=(z>0)?0:18;
 885         xyz[0]=x;
 886         xyz[1]=y;
 887         xyz[2]=z;
 888         Double_t phi  = TMath::ATan2(y,x);
 889         DistortPoint(xyz,roc);
 890         Double_t r1    = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 891         Double_t phi1  = TMath::ATan2(xyz[1],xyz[0]);
 892         if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
 893         if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
 894         Double_t dx = xyz[0]-x;
 895         Double_t dy = xyz[1]-y;
 896         Double_t dz = xyz[2]-z;
 897         Double_t dr=r1-r;
 898         Double_t drphi=(phi1-phi)*r;
 899         (*pcstream)<<"distortion"<<
 900           "x="<<x<<           // original position
 901           "y="<<y<<
 902           "z="<<z<<
 903           "r="<<r<<
 904           "phi="<<phi<<
 905           "x1="<<xyz[0]<<      // distorted position
 906           "y1="<<xyz[1]<<
 907           "z1="<<xyz[2]<<
 908           "r1="<<r1<<
 909           "phi1="<<phi1<<
 910           //
 911           "dx="<<dx<<          // delta position
 912           "dy="<<dy<<
 913           "dz="<<dz<<
 914           "dr="<<dr<<
 915           "drphi="<<drphi<<
 916           "\n";
 917       }
 918     }
 919   }
 920   delete pcstream;
 921   TFile f(Form("correction%s.root",GetName()));
 922   TTree * tree = (TTree*)f.Get("distortion");
 923   TTree * tree2= tree->CopyTree("1");
 924   tree2->SetName(Form("dist%s",GetName()));
 925   tree2->SetDirectory(0);
 926   delete tree;
 927   return tree2;
 928 }
 929
 930
 931
 932
 933 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
 934   //
 935   // Make a fit tree:
 936   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
 937   // calculates partial distortions
 938   // Partial distortion is stored in the resulting tree
 939   // Output is storred in the file distortion_<dettype>_<partype>.root
 940   // Partial  distortion is stored with the name given by correction name
 941   //
 942   //
 943   // Parameters of function:
 944   // input     - input tree
 945   // dtype     - distortion type 0 - ITSTPC,  1 -TPCTRD, 2 - TPCvertex
 946   // ppype     - parameter type
 947   // corrArray - array with partial corrections
 948   // step      - skipe entries  - if 1 all entries processed - it is slow
 949   // debug     0 if debug on also space points dumped - it is slow
 950   const Double_t kMaxSnp = 0.85;
 951   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 952   //  const Double_t kB2C=-0.299792458e-3;
 953   const Int_t kMinEntries=50;
 954   Double_t phi,theta, snp, mean,rms, entries;
 955   tinput->SetBranchAddress("theta",&theta);
 956   tinput->SetBranchAddress("phi", &phi);
 957   tinput->SetBranchAddress("snp",&snp);
 958   tinput->SetBranchAddress("mean",&mean);
 959   tinput->SetBranchAddress("rms",&rms);
 960   tinput->SetBranchAddress("entries",&entries);
 961   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
 962   //
 963   Int_t nentries=tinput->GetEntries();
 964   Int_t ncorr=corrArray->GetEntries();
 965   Double_t corrections[100]={0}; //
 966   Double_t tPar[5];
 967   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 968   Double_t refX=0;
 969   Int_t dir=0;
 970   if (dtype==0) {refX=85.; dir=-1;}
 971   if (dtype==1) {refX=275.; dir=1;}
 972   if (dtype==2) {refX=85.; dir=-1;}
 973   if (dtype==3) {refX=360.; dir=-1;}
 974   //
 975   for (Int_t ientry=0; ientry<nentries; ientry+=step){
 976     tinput->GetEntry(ientry);
 977     if (TMath::Abs(snp)>kMaxSnp) continue;
 978     tPar[0]=0;
 979     tPar[1]=theta*refX;
 980     tPar[2]=snp;
 981     tPar[3]=theta;
 982     tPar[4]=(gRandom->Rndm()-0.5)*0.02;  // should be calculated - non equal to 0
 983     Double_t bz=AliTrackerBase::GetBz();
 984     if (refX>10. && TMath::Abs(bz)>0.1 )  tPar[4]=snp/(refX*bz*kB2C*2);
 985     tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
 986     AliExternalTrackParam track(refX,phi,tPar,cov);
 987     Double_t xyz[3];
 988     track.GetXYZ(xyz);
 989     Int_t id=0;
 990     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
 991     if (ptype==4 &&bz<0) mean*=-1;  // interpret as curvature
 992     (*pcstream)<<"fit"<<
 993       "bz="<<bz<<         // magnetic filed used
 994       "dtype="<<dtype<<   // detector match type
 995       "ptype="<<ptype<<   // parameter type
 996       "theta="<<theta<<   // theta
 997       "phi="<<phi<<       // phi
 998       "snp="<<snp<<       // snp
 999       "mean="<<mean<<     // mean dist value
1000       "rms="<<rms<<       // rms
1001       "gx="<<xyz[0]<<         // global position at reference
1002       "gy="<<xyz[1]<<         // global position at reference
1003       "gz="<<xyz[2]<<         // global position at reference
1004       "dRrec="<<dRrec<<      // delta Radius in reconstruction
1005       "id="<<id<<             // track id
1006       "entries="<<entries;// number of entries in bin
1007     //
1008     for (Int_t icorr=0; icorr<ncorr; icorr++) {
1009       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1010       corrections[icorr]=0;
1011       if (entries>kMinEntries){
1012         AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1013         AliExternalTrackParam *trackOut = 0;
1014         if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1015         if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1016         if (dtype==0) {refX=85.; dir=-1;}
1017         if (dtype==1) {refX=275.; dir=1;}
1018         if (dtype==2) {refX=0; dir=-1;}
1019         if (dtype==3) {refX=360.; dir=-1;}
1020         //
1021         if (trackOut){
1022           AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1023           trackOut->Rotate(trackIn.GetAlpha());
1024           trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1025           //
1026           corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1027           delete trackOut;
1028         }else{
1029           corrections[icorr]=0;
1030         }
1031         if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature
1032       }
1033       Double_t dRdummy=0;
1034       (*pcstream)<<"fit"<<
1035         Form("%s=",corr->GetName())<<corrections[icorr]<<   // dump correction value
1036         Form("dR%s=",corr->GetName())<<dRdummy;   // dump dummy correction value not needed for tracks
1037                                                   // for points it is neccessary
1038     }
1039     (*pcstream)<<"fit"<<"\n";
1040   }
1041   delete pcstream;
1042 }
1043
1044
1045
1046 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1047   //
1048   // Make a laser fit tree for global minimization
1049   //
1050   const Double_t cutErrY=0.1;
1051   const Double_t cutErrZ=0.1;
1052   const Double_t kEpsilon=0.00000001;
1053   TVectorD *vecdY=0;
1054   TVectorD *vecdZ=0;
1055   TVectorD *veceY=0;
1056   TVectorD *veceZ=0;
1057   AliTPCLaserTrack *ltr=0;
1058   AliTPCLaserTrack::LoadTracks();
1059   tree->SetBranchAddress("dY.",&vecdY);
1060   tree->SetBranchAddress("dZ.",&vecdZ);
1061   tree->SetBranchAddress("eY.",&veceY);
1062   tree->SetBranchAddress("eZ.",&veceZ);
1063   tree->SetBranchAddress("LTr.",&ltr);
1064   Int_t entries= tree->GetEntries();
1065   TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1066   Double_t bz=AliTrackerBase::GetBz();
1067   //
1068
1069   for (Int_t ientry=0; ientry<entries; ientry++){
1070     tree->GetEntry(ientry);
1071     if (!ltr->GetVecGX()){
1072       ltr->UpdatePoints();
1073     }
1074     TVectorD * delta= (itype==0)? vecdY:vecdZ;
1075     TVectorD * err= (itype==0)? veceY:veceZ;
1076
1077     for (Int_t irow=0; irow<159; irow++){
1078       Int_t nentries = 1000;
1079       if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1080       if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1081       Int_t dtype=4;
1082       Double_t phi   =(*ltr->GetVecPhi())[irow];
1083       Double_t theta =ltr->GetTgl();
1084       Double_t mean=delta->GetMatrixArray()[irow];
1085       Double_t gx=0,gy=0,gz=0;
1086       Double_t snp = (*ltr->GetVecP2())[irow];
1087       Double_t rms = 0.1+err->GetMatrixArray()[irow];
1088       gx = (*ltr->GetVecGX())[irow];
1089       gy = (*ltr->GetVecGY())[irow];
1090       gz = (*ltr->GetVecGZ())[irow];
1091       Int_t bundle= ltr->GetBundle();
1092       Double_t dRrec=0;
1093       //
1094       // get delta R used in reconstruction
1095       AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
1096       AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1097       const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1098       Double_t xyz0[3]={gx,gy,gz};
1099       Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1100       //
1101       // old ExB correction
1102       //
1103       if(recoParam&&recoParam->GetUseExBCorrection()) {
1104         Double_t xyz1[3]={gx,gy,gz};
1105         calib->GetExB()->Correct(xyz0,xyz1);
1106         Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1107         dRrec=oldR-newR;
1108       }
1109       if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1110         Float_t xyz1[3]={gx,gy,gz};
1111         Int_t sector=(gz>0)?0:18;
1112         correction->CorrectPoint(xyz1, sector);
1113         Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1114         dRrec=oldR-newR;
1115       }
1116
1117
1118       (*pcstream)<<"fit"<<
1119         "bz="<<bz<<         // magnetic filed used
1120         "dtype="<<dtype<<   // detector match type
1121         "ptype="<<itype<<   // parameter type
1122         "theta="<<theta<<   // theta
1123         "phi="<<phi<<       // phi
1124         "snp="<<snp<<       // snp
1125         "mean="<<mean<<     // mean dist value
1126         "rms="<<rms<<       // rms
1127         "gx="<<gx<<         // global position
1128         "gy="<<gy<<         // global position
1129         "gz="<<gz<<         // global position
1130         "dRrec="<<dRrec<<      // delta Radius in reconstruction
1131         "id="<<bundle<<     //bundle
1132         "entries="<<nentries;// number of entries in bin
1133       //
1134       //
1135       Double_t ky = TMath::Tan(TMath::ASin(snp));
1136       Int_t ncorr = corrArray->GetEntries();
1137       Double_t r0   = TMath::Sqrt(gx*gx+gy*gy);
1138       Double_t phi0 = TMath::ATan2(gy,gx);
1139       Double_t distortions[1000]={0};
1140       Double_t distortionsR[1000]={0};
1141       for (Int_t icorr=0; icorr<ncorr; icorr++) {
1142         AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1143         Float_t distPoint[3]={gx,gy,gz};
1144         Int_t sector= (gz>0)? 0:18;
1145         if (r0>80){
1146           corr->DistortPoint(distPoint, sector);
1147         }
1148         // Double_t value=distPoint[2]-gz;
1149         if (itype==0){
1150           Double_t r1   = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1151           Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1152           Double_t drphi= r0*(phi1-phi0);
1153           Double_t dr   = r1-r0;
1154           distortions[icorr]  = drphi-ky*dr;
1155           distortionsR[icorr] = dr;
1156         }
1157         (*pcstream)<<"fit"<<
1158           Form("%s=",corr->GetName())<<distortions[icorr]<<    // dump correction value
1159           Form("dR%s=",corr->GetName())<<distortionsR[icorr];   // dump correction R  value
1160       }
1161       (*pcstream)<<"fit"<<"\n";
1162     }
1163   }
1164   delete pcstream;
1165 }
1166
1167
1168
1169 void   AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
1170   //
1171   // make a distortion map out ou fthe residual histogram
1172   // Results are written to the debug streamer - pcstream
1173   // Parameters:
1174   //   his0       - input (4D) residual histogram
1175   //   pcstream   - file to write the tree
1176   //   run        - run number
1177   // marian.ivanov@cern.ch
1178   const Int_t kMinEntries=50;
1179   Int_t nbins1=his0->GetAxis(1)->GetNbins();
1180   Int_t first1=his0->GetAxis(1)->GetFirst();
1181   Int_t last1 =his0->GetAxis(1)->GetLast();
1182   //
1183   Double_t bz=AliTrackerBase::GetBz();
1184   Int_t idim[4]={0,1,2,3};
1185   for (Int_t ibin1=first1; ibin1<last1; ibin1++){   //axis 1 - theta
1186     //
1187     his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1188     Double_t       x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1189     THnSparse * his1 = his0->Projection(4,idim);  // projected histogram according range1
1190     Int_t nbins3     = his1->GetAxis(3)->GetNbins();
1191     Int_t first3     = his1->GetAxis(3)->GetFirst();
1192     Int_t last3      = his1->GetAxis(3)->GetLast();
1193     //
1194
1195     for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){   // axis 3 - local angle
1196       his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1197       Double_t      x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1198       if (ibin3<first3) {
1199         his1->GetAxis(3)->SetRangeUser(-1,1);
1200         x3=0;
1201       }
1202       THnSparse * his3= his1->Projection(4,idim);         //projected histogram according selection 3
1203       Int_t nbins2    = his3->GetAxis(2)->GetNbins();
1204       Int_t first2    = his3->GetAxis(2)->GetFirst();
1205       Int_t last2     = his3->GetAxis(2)->GetLast();
1206       //
1207       for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1208         his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1209         Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1210         TH1 * hisDelta = his3->Projection(0);
1211         //
1212         Double_t entries = hisDelta->GetEntries();
1213         Double_t mean=0, rms=0;
1214         if (entries>kMinEntries){
1215           mean    = hisDelta->GetMean();
1216           rms = hisDelta->GetRMS();
1217         }
1218         (*pcstream)<<hname<<
1219           "run="<<run<<
1220           "bz="<<bz<<
1221           "theta="<<x1<<
1222           "phi="<<x2<<
1223           "snp="<<x3<<
1224           "entries="<<entries<<
1225           "mean="<<mean<<
1226           "rms="<<rms<<
1227           "\n";
1228         delete hisDelta;
1229         printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1230       }
1231       delete his3;
1232     }
1233     delete his1;
1234   }
1235 }
1236
1237
1238
1239
1240
1241 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1242   //
1243   // Store object in the OCDB
1244   // By default the object is stored in the current directory
1245   // default comment consit of user name and the date
1246   //
1247   TString ocdbStorage="";
1248   ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1249   AliCDBMetaData *metaData= new AliCDBMetaData();
1250   metaData->SetObjectClassName("AliTPCCorrection");
1251   metaData->SetResponsible("Marian Ivanov");
1252   metaData->SetBeamPeriod(1);
1253   metaData->SetAliRootVersion("05-25-01"); //root version
1254   TString userName=gSystem->GetFromPipe("echo $USER");
1255   TString date=gSystem->GetFromPipe("date");
1256
1257   if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1258   if (comment) metaData->SetComment(comment);
1259   AliCDBId* id1=NULL;
1260   id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1261   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1262   gStorage->Put(this, (*id1), metaData);
1263 }
1264
1265
1266 void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
1267
1268   AliVertexerTracks *vertexer = new AliVertexerTracks(5);// 5kGaus
1269
1270   TObjArray   ATrk;               // Original Track array of Aside
1271   TObjArray   dATrk;              // Distorted Track array of A side
1272   UShort_t    *AId = new UShort_t[nTracks];      // A side Track ID
1273   TObjArray   CTrk;
1274   TObjArray   dCTrk;
1275   UShort_t    *CId = new UShort_t [nTracks];
1276   Int_t ID=0;
1277   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1278   TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
1279   fpt.SetParameters(7.24,0.120);
1280   fpt.SetNpx(10000);
1281   for(Int_t nt=0; nt<nTracks; nt++){
1282     Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
1283     Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
1284     Double_t pt = fpt.GetRandom();// momentum for f1
1285     Short_t sign=1;
1286     if(gRandom->Rndm() < 0.5){
1287       sign =1;
1288     }else{
1289       sign=-1;
1290     }
1291
1292     Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1293     Double_t pxyz[3];
1294     pxyz[0]=pt*TMath::Cos(phi);
1295     pxyz[1]=pt*TMath::Sin(phi);
1296     pxyz[2]=pt*TMath::Tan(theta);
1297     Double_t cv[21]={0};
1298     AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
1299
1300     Double_t refX=1.;
1301     Int_t dir=-1;
1302     AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir,  NULL);
1303     if (!td) continue;
1304     if (pcstream) (*pcstream)<<"track"<<
1305       "eta="<<eta<<
1306       "theta="<<theta<<
1307       "tOrig.="<<t<<
1308       "td.="<<td<<
1309       "\n";
1310     if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
1311       if (td){
1312         dATrk.AddLast(td);
1313         ATrk.AddLast(t);
1314         Int_t nn=ATrk.GetEntriesFast();
1315         AId[nn]=ID;
1316       }
1317     }else if(( eta<-0.07 )&&( eta>-etaCuts )){
1318       if (td){
1319         dCTrk.AddLast(td);
1320         CTrk.AddLast(t);
1321         Int_t nn=CTrk.GetEntriesFast();
1322         CId[nn]=ID;
1323       }
1324     }
1325     ID++;
1326   }// end of track loop
1327
1328   vertexer->SetTPCMode();
1329   vertexer->SetConstraintOff();
1330
1331   aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dATrk,AId));
1332   avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&ATrk,AId));
1333   cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dCTrk,CId));
1334   cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&CTrk,CId));
1335   if (pcstream) (*pcstream)<<"vertex"<<
1336     "x="<<orgVertex[0]<<
1337     "y="<<orgVertex[1]<<
1338     "z="<<orgVertex[2]<<
1339     "av.="<<&aV<<              // distorted vertex A side
1340     "cv.="<<&cV<<              // distroted vertex C side
1341     "avO.="<<&avOrg<<         // original vertex A side
1342     "cvO.="<<&cvOrg<<
1343     "\n";
1344   delete []AId;
1345   delete []CId;
1346 }