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
  42 #include <TH2F.h>
  43 #include <TMath.h>
  44 #include <TROOT.h>
  45 #include <TTreeStream.h>
  46 #include <TTimeStamp.h>
  47
  48 #include  "AliExternalTrackParam.h"
  49 #include  "AliTrackPointArray.h"
  50 #include  "TDatabasePDG.h"
  51 #include  "AliTrackerBase.h"
  52 #include  "AliTPCROC.h"
  53
  54
  55 #include "AliTPCCorrection.h"
  56
  57 ClassImp(AliTPCCorrection)
  58
  59 // FIXME: the following values should come from the database
  60 const Double_t AliTPCCorrection::fgkTPC_Z0   =249.7;     // nominal gating grid position
  61 const Double_t AliTPCCorrection::fgkIFCRadius= 83.06;    // Mean Radius of the Inner Field Cage ( 82.43 min,  83.70 max) (cm)
  62 const Double_t AliTPCCorrection::fgkOFCRadius=254.5;     // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
  63 const Double_t AliTPCCorrection::fgkZOffSet  = 0.2;      // Offset from CE: calculate all distortions closer to CE as if at this point
  64 const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
  65 const Double_t AliTPCCorrection::fgkGG       =-70.0;     // Gating Grid voltage (volts)
  66
  67
  68 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
  69 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] =  {
  70 84.0,   84.5,   85.0,   85.5,   86.0,  87.0,    88.0,
  71 90.0,   92.0,   94.0,   96.0,   98.0,  100.0,  102.0,  104.0,  106.0,  108.0,
  72 110.0,  112.0,  114.0,  116.0,  118.0,  120.0,  122.0,  124.0,  126.0,  128.0,
  73 130.0,  132.0,  134.0,  136.0,  138.0,  140.0,  142.0,  144.0,  146.0,  148.0,
  74 150.0,  152.0,  154.0,  156.0,  158.0,  160.0,  162.0,  164.0,  166.0,  168.0,
  75 170.0,  172.0,  174.0,  176.0,  178.0,  180.0,  182.0,  184.0,  186.0,  188.0,
  76 190.0,  192.0,  194.0,  196.0,  198.0,  200.0,  202.0,  204.0,  206.0,  208.0,
  77 210.0,  212.0,  214.0,  216.0,  218.0,  220.0,  222.0,  224.0,  226.0,  228.0,
  78 230.0,  232.0,  234.0,  236.0,  238.0,  240.0,  242.0,  244.0,  246.0,  248.0,
  79 249.0,  249.5,  250.0,  251.5,  252.0  } ;
  80
  81 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ]     =   {
  82 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
  83 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
  84 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
  85 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
  86 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
  87 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
  88 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
  89 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
  90 -100.0,  -98.0,  -96.0,  -94.0,  -92.0,  -90.0,  -88.0,  -86.0,  -84.0,  -82.0,
  91 -80.0,  -78.0,  -76.0,  -74.0,  -72.0,  -70.0,  -68.0,  -66.0,  -64.0,  -62.0,
  92 -60.0,  -58.0,  -56.0,  -54.0,  -52.0,  -50.0,  -48.0,  -46.0,  -44.0,  -42.0,
  93 -40.0,  -38.0,  -36.0,  -34.0,  -32.0,  -30.0,  -28.0,  -26.0,  -24.0,  -22.0,
  94 -20.0,  -18.0,  -16.0,  -14.0,  -12.0,  -10.0,   -8.0,   -6.0,   -4.0,   -2.0,
  95 -1.0,   -0.5,   -0.2,   -0.1,  -0.05,   0.05,    0.1,    0.2,    0.5,    1.0,
  96  2.0,    4.0,    6.0,    8.0,   10.0,   12.0,   14.0,   16.0,   18.0,   20.0,
  97  22.0,   24.0,   26.0,   28.0,   30.0,   32.0,   34.0,   36.0,   38.0,   40.0,
  98  42.0,   44.0,   46.0,   48.0,   50.0,   52.0,   54.0,   56.0,   58.0,   60.0,
  99  62.0,   64.0,   66.0,   68.0,   70.0,   72.0,   74.0,   76.0,   78.0,   80.0,
 100  82.0,   84.0,   86.0,   88.0,   90.0,   92.0,   94.0,   96.0,   98.0,  100.0,
 101 102.0,  104.0,  106.0,  108.0,  110.0,  112.0,  114.0,  116.0,  118.0,  120.0,
 102 122.0,  124.0,  126.0,  128.0,  130.0,  132.0,  134.0,  136.0,  138.0,  140.0,
 103 142.0,  144.0,  146.0,  148.0,  150.0,  152.0,  154.0,  156.0,  158.0,  160.0,
 104 162.0,  164.0,  166.0,  168.0,  170.0,  172.0,  174.0,  176.0,  178.0,  180.0,
 105 182.0,  184.0,  186.0,  188.0,  190.0,  192.0,  194.0,  196.0,  198.0,  200.0,
 106 202.0,  204.0,  206.0,  208.0,  210.0,  212.0,  214.0,  216.0,  218.0,  220.0,
 107 222.0,  224.0,  226.0,  228.0,  230.0,  232.0,  234.0,  236.0,  238.0,  240.0,
 108 242.0,  243.0,  244.0,  245.0,  246.0,  247.0,  248.0,  248.5,  249.0,  249.5   } ;
 109
 110
 111
 112 AliTPCCorrection::AliTPCCorrection()
 113   : TNamed("correction_unity","unity"),fJLow(0),fKLow(0)
 114 {
 115   //
 116   // default constructor
 117   //
 118 }
 119
 120 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
 121   : TNamed(name,title),fJLow(0),fKLow(0)
 122 {
 123   //
 124   // default constructor, that set the name and title
 125   //
 126 }
 127
 128 AliTPCCorrection::~AliTPCCorrection() {
 129   //
 130   // virtual destructor
 131   //
 132 }
 133
 134 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
 135   //
 136   // Corrects the initial coordinates x (cartesian coordinates)
 137   // according to the given effect (inherited classes)
 138   // roc represents the TPC read out chamber (offline numbering convention)
 139   //
 140   Float_t dx[3];
 141   GetCorrection(x,roc,dx);
 142   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 143 }
 144
 145 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 146   //
 147   // Corrects the initial coordinates x (cartesian coordinates) and stores the new
 148   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 149   // roc represents the TPC read out chamber (offline numbering convention)
 150   //
 151   Float_t dx[3];
 152   GetCorrection(x,roc,dx);
 153   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 154 }
 155
 156 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
 157   //
 158   // Distorts the initial coordinates x (cartesian coordinates)
 159   // according to the given effect (inherited classes)
 160   // roc represents the TPC read out chamber (offline numbering convention)
 161   //
 162   Float_t dx[3];
 163   GetDistortion(x,roc,dx);
 164   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 165 }
 166
 167 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 168   //
 169   // Distorts the initial coordinates x (cartesian coordinates) and stores the new
 170   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 171   // roc represents the TPC read out chamber (offline numbering convention)
 172   //
 173   Float_t dx[3];
 174   GetDistortion(x,roc,dx);
 175   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 176 }
 177
 178 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
 179   //
 180   // This function delivers the correction values dx in respect to the inital coordinates x
 181   // roc represents the TPC read out chamber (offline numbering convention)
 182   // Note: The dx is overwritten by the inherited effectice class ...
 183   //
 184   for (Int_t j=0;j<3;++j) { dx[j]=0.; }
 185 }
 186
 187 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
 188   //
 189   // This function delivers the distortion values dx in respect to the inital coordinates x
 190   // roc represents the TPC read out chamber (offline numbering convention)
 191   //
 192   GetCorrection(x,roc,dx);
 193   for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
 194 }
 195
 196 void AliTPCCorrection::Init() {
 197   //
 198   // Initialization funtion (not used at the moment)
 199   //
 200 }
 201
 202 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
 203   //
 204   // Update function
 205   //
 206 }
 207
 208 void AliTPCCorrection::Print(Option_t* /*option*/) const {
 209   //
 210   // Print function to check which correction classes are used
 211   // option=="d" prints details regarding the setted magnitude
 212   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 213   //
 214   printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
 215 }
 216
 217 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t /*t1*/,Float_t /*t2*/) {
 218   //
 219   // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
 220   // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
 221   // calibration run
 222   //
 223   // SetOmegaTauT1T2(omegaTau, t1, t2);
 224 }
 225
 226 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
 227   //
 228   // Simple plot functionality.
 229   // Returns a 2d hisogram which represents the corrections in radial direction (dr)
 230   // in respect to position z within the XY plane.
 231   // The histogramm has nx times ny entries.
 232   //
 233
 234   TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
 235                      nx,-250.,250.,ny,-250.,250.);
 236   Float_t x[3],dx[3];
 237   x[2]=z;
 238   Int_t roc=z>0.?0:18; // FIXME
 239   for (Int_t iy=1;iy<=ny;++iy) {
 240     x[1]=h->GetYaxis()->GetBinCenter(iy);
 241     for (Int_t ix=1;ix<=nx;++ix) {
 242       x[0]=h->GetXaxis()->GetBinCenter(ix);
 243       GetCorrection(x,roc,dx);
 244       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 245       if (90.<=r0 && r0<=250.) {
 246         Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 247         h->SetBinContent(ix,iy,r1-r0);
 248       }
 249       else
 250         h->SetBinContent(ix,iy,0.);
 251     }
 252   }
 253   return h;
 254 }
 255
 256 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
 257   //
 258   // Simple plot functionality.
 259   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 260   // in respect to position z within the XY plane.
 261   // The histogramm has nx times ny entries.
 262   //
 263
 264   TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
 265                      nx,-250.,250.,ny,-250.,250.);
 266   Float_t x[3],dx[3];
 267   x[2]=z;
 268   Int_t roc=z>0.?0:18; // FIXME
 269   for (Int_t iy=1;iy<=ny;++iy) {
 270     x[1]=h->GetYaxis()->GetBinCenter(iy);
 271     for (Int_t ix=1;ix<=nx;++ix) {
 272       x[0]=h->GetXaxis()->GetBinCenter(ix);
 273       GetCorrection(x,roc,dx);
 274       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 275       if (90.<=r0 && r0<=250.) {
 276         Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 277         Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 278
 279         Float_t dphi=phi1-phi0;
 280         if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 281         if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 282
 283         h->SetBinContent(ix,iy,r0*dphi);
 284       }
 285       else
 286         h->SetBinContent(ix,iy,0.);
 287     }
 288   }
 289   return h;
 290 }
 291
 292 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
 293   //
 294   // Simple plot functionality.
 295   // Returns a 2d hisogram which represents the corrections in r direction (dr)
 296   // in respect to angle phi within the ZR plane.
 297   // The histogramm has nx times ny entries.
 298   //
 299   TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
 300                      nz,-250.,250.,nr,85.,250.);
 301   Float_t x[3],dx[3];
 302   for (Int_t ir=1;ir<=nr;++ir) {
 303     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 304     x[0]=radius*TMath::Cos(phi);
 305     x[1]=radius*TMath::Sin(phi);
 306     for (Int_t iz=1;iz<=nz;++iz) {
 307       x[2]=h->GetXaxis()->GetBinCenter(iz);
 308       Int_t roc=x[2]>0.?0:18; // FIXME
 309       GetCorrection(x,roc,dx);
 310       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 311       Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 312       h->SetBinContent(iz,ir,r1-r0);
 313     }
 314   }
 315   printf("SDF\n");
 316   return h;
 317
 318 }
 319
 320 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
 321   //
 322   // Simple plot functionality.
 323   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 324   // in respect to angle phi within the ZR plane.
 325   // The histogramm has nx times ny entries.
 326   //
 327   TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
 328                      nz,-250.,250.,nr,85.,250.);
 329   Float_t x[3],dx[3];
 330   for (Int_t iz=1;iz<=nz;++iz) {
 331     x[2]=h->GetXaxis()->GetBinCenter(iz);
 332     Int_t roc=x[2]>0.?0:18; // FIXME
 333     for (Int_t ir=1;ir<=nr;++ir) {
 334       Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 335       x[0]=radius*TMath::Cos(phi);
 336       x[1]=radius*TMath::Sin(phi);
 337       GetCorrection(x,roc,dx);
 338       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 339       Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 340       Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 341
 342       Float_t dphi=phi1-phi0;
 343       if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 344       if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 345
 346       h->SetBinContent(iz,ir,r0*dphi);
 347     }
 348   }
 349   return h;
 350 }
 351
 352 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
 353                                    const char *xlabel,const char *ylabel,const char *zlabel,
 354                                   Int_t nbinsx,Double_t xlow,Double_t xup,
 355                                   Int_t nbinsy,Double_t ylow,Double_t yup) {
 356   //
 357   // Helper function to create a 2d histogramm of given size
 358   //
 359
 360   TString hname=name;
 361   Int_t i=0;
 362   if (gDirectory) {
 363     while (gDirectory->FindObject(hname.Data())) {
 364       hname =name;
 365       hname+="_";
 366       hname+=i;
 367       ++i;
 368     }
 369   }
 370   TH2F *h=new TH2F(hname.Data(),title,
 371                    nbinsx,xlow,xup,
 372                    nbinsy,ylow,yup);
 373   h->GetXaxis()->SetTitle(xlabel);
 374   h->GetYaxis()->SetTitle(ylabel);
 375   h->GetZaxis()->SetTitle(zlabel);
 376   h->SetStats(0);
 377   return h;
 378 }
 379
 380
 381 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
 382
 383 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
 384                                                   const Double_t er[kNZ][kNR], Double_t &er_value )
 385 {
 386   //
 387   // Interpolate table - 2D interpolation
 388   //
 389   Double_t save_er[10] ;
 390
 391   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
 392   Search( kNR,   fgkRList,  r,   fKLow   ) ;
 393   if ( fJLow < 0 ) fJLow = 0 ;   // check if out of range
 394   if ( fKLow < 0 ) fKLow = 0 ;
 395   if ( fJLow + order  >=    kNZ - 1 ) fJLow =   kNZ - 1 - order ;
 396   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 397
 398   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
 399       save_er[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
 400   }
 401   er_value = Interpolate( &fgkZList[fJLow], save_er, order, z )   ;
 402
 403 }
 404
 405
 406 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
 407                                        const Int_t order, const Double_t x )
 408 {
 409   //
 410   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
 411   //
 412
 413   Double_t y ;
 414   if ( order == 2 ) {                // Quadratic Interpolation = 2
 415     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
 416     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
 417     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
 418   } else {                           // Linear Interpolation = 1
 419     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
 420   }
 421
 422   return (y);
 423
 424 }
 425
 426
 427 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low )
 428 {
 429   //
 430   // Search an ordered table by starting at the most recently used point
 431   //
 432
 433   Long_t middle, high ;
 434   Int_t  ascend = 0, increment = 1 ;
 435
 436   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
 437
 438   if ( low < 0 || low > n-1 ) {
 439     low = -1 ; high = n ;
 440   } else {                                            // Ordered Search phase
 441     if ( (Int_t)( x >= xArray[low] ) == ascend )  {
 442       if ( low == n-1 ) return ;
 443       high = low + 1 ;
 444       while ( (Int_t)( x >= xArray[high] ) == ascend ) {
 445         low = high ;
 446         increment *= 2 ;
 447         high = low + increment ;
 448         if ( high > n-1 )  {  high = n ; break ;  }
 449       }
 450     } else {
 451       if ( low == 0 )  {  low = -1 ;  return ;  }
 452       high = low - 1 ;
 453       while ( (Int_t)( x < xArray[low] ) == ascend ) {
 454         high = low ;
 455         increment *= 2 ;
 456         if ( increment >= high )  {  low = -1 ;  break ;  }
 457         else  low = high - increment ;
 458       }
 459     }
 460   }
 461
 462   while ( (high-low) != 1 ) {                     // Binary Search Phase
 463     middle = ( high + low ) / 2 ;
 464     if ( (Int_t)( x >= xArray[middle] ) == ascend )
 465       low = middle ;
 466     else
 467       high = middle ;
 468   }
 469
 470   if ( x == xArray[n-1] ) low = n-2 ;
 471   if ( x == xArray[0]   ) low = 0 ;
 472
 473 }
 474
 475
 476 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(const AliExternalTrackParam * trackIn, Double_t refX, Int_t dir,TTreeSRedirector *pcstream){
 477   //
 478   // Fit the track parameters - without and with distortion
 479   // 1. Space points in the TPC are simulated along the trajectory
 480   // 2. Space points distorted
 481   // 3. Fits  the non distorted and distroted track to the reference plane at refX
 482   // 4. For visualization and debugging  purposes the space points and tracks can be stored  in the tree - using the TTreeSRedirector functionality
 483   //
 484   // trackIn   - input track parameters
 485   // refX     - reference X to fit the track
 486   // dir      - direction - out=1 or in=-1
 487   // pcstream -  debug streamer to check the results
 488   //
 489   // see AliExternalTrackParam.h documentation:
 490   // track1.fP[0] - local y (rphi)
 491   // track1.fP[1] - z
 492   // track1.fP[2] - sinus of local inclination angle
 493   // track1.fP[3] - tangent of deep angle
 494   // track1.fP[4] - 1/pt
 495   AliTPCROC * roc = AliTPCROC::Instance();
 496   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
 497   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
 498   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 499
 500   const Double_t kMaxSnp = 0.85;
 501   const Double_t kSigmaY=0.1;
 502   const Double_t kSigmaZ=0.1;
 503   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 504
 505   AliExternalTrackParam  track(*trackIn); //
 506   // generate points
 507   AliTrackPointArray pointArray0(npoints0);
 508   AliTrackPointArray pointArray1(npoints0);
 509   Double_t xyz[3];
 510   AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
 511   //
 512   // simulate the track
 513   Int_t npoints=0;
 514   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
 515   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
 516     AliTrackerBase::PropagateTrackTo(&track,radius,kMass,3,kTRUE,kMaxSnp);
 517     track.GetXYZ(xyz);
 518     AliTrackPoint pIn0;                               // space point
 519     AliTrackPoint pIn1;
 520     Int_t roc= (xyz[2]>0)? 0:18;
 521     pointArray0.GetPoint(pIn0,npoints);
 522     pointArray1.GetPoint(pIn1,npoints);
 523     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 524     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
 525     DistortPoint(distPoint, roc);
 526     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
 527     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
 528     //
 529     track.Rotate(alpha);
 530     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
 531     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
 532     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
 533     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
 534     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
 535     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
 536     pointArray0.AddPoint(npoints, &pIn0);
 537     pointArray1.AddPoint(npoints, &pIn1);
 538     npoints++;
 539     if (npoints>=npoints0) break;
 540   }
 541   //
 542   // refit track
 543   //
 544   AliExternalTrackParam *track0=0;
 545   AliExternalTrackParam *track1=0;
 546   AliTrackPoint   point1,point2,point3;
 547   if (dir==1) {  //make seed inner
 548     pointArray0.GetPoint(point1,1);
 549     pointArray0.GetPoint(point2,10);
 550     pointArray0.GetPoint(point3,20);
 551   }
 552   if (dir==-1){ //make seed outer
 553     pointArray0.GetPoint(point1,npoints-20);
 554     pointArray0.GetPoint(point2,npoints-10);
 555     pointArray0.GetPoint(point3,npoints-1);
 556   }
 557   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
 558   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
 559
 560
 561   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
 562     Int_t ipoint= (dir>0) ? ipoint: npoints-1-jpoint;
 563     //
 564     AliTrackPoint pIn0;
 565     AliTrackPoint pIn1;
 566     pointArray0.GetPoint(pIn0,ipoint);
 567     pointArray1.GetPoint(pIn1,ipoint);
 568     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
 569     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
 570     //
 571     AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,1,kFALSE,kMaxSnp);
 572     AliTrackerBase::PropagateTrackTo(track1,prot1.GetX(),kMass,1,kFALSE,kMaxSnp);
 573     track.GetXYZ(xyz);
 574     //
 575     Double_t pointPos[2]={0,0};
 576     Double_t pointCov[3]={0,0,0};
 577     pointPos[0]=prot0.GetY();//local y
 578     pointPos[1]=prot0.GetZ();//local z
 579     pointCov[0]=prot0.GetCov()[3];//simay^2
 580     pointCov[1]=prot0.GetCov()[4];//sigmayz
 581     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
 582     track0->Update(pointPos,pointCov);
 583     //
 584     pointPos[0]=prot1.GetY();//local y
 585     pointPos[1]=prot1.GetZ();//local z
 586     pointCov[0]=prot1.GetCov()[3];//simay^2
 587     pointCov[1]=prot1.GetCov()[4];//sigmayz
 588     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
 589     track1->Update(pointPos,pointCov);
 590   }
 591
 592   AliTrackerBase::PropagateTrackTo(track0,refX,kMass,2.,kTRUE,kMaxSnp);
 593   track1->Rotate(track0->GetAlpha());
 594   AliTrackerBase::PropagateTrackTo(track1,refX,kMass,2.,kFALSE,kMaxSnp);
 595
 596   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
 597     "point0.="<<&pointArray0<<   //  points
 598     "point1.="<<&pointArray1<<   //  distorted points
 599     "trackIn.="<<&track<<       //  original track
 600     "track0.="<<track0<<         //  fitted track
 601     "track1.="<<track1<<         //  fitted distorted track
 602     "\n";
 603   delete track0;
 604   return track1;
 605 }
 606
 607