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
  54 #include  "AliExternalTrackParam.h"
  55 #include  "AliTrackPointArray.h"
  56 #include  "TDatabasePDG.h"
  57 #include  "AliTrackerBase.h"
  58 #include  "AliTPCROC.h"
  59
  60
  61 #include "AliTPCCorrection.h"
  62
  63 ClassImp(AliTPCCorrection)
  64
  65 // FIXME: the following values should come from the database
  66 const Double_t AliTPCCorrection::fgkTPC_Z0   =249.7;     // nominal gating grid position
  67 const Double_t AliTPCCorrection::fgkIFCRadius= 83.06;    // Mean Radius of the Inner Field Cage ( 82.43 min,  83.70 max) (cm)
  68 const Double_t AliTPCCorrection::fgkOFCRadius=254.5;     // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
  69 const Double_t AliTPCCorrection::fgkZOffSet  = 0.2;      // Offset from CE: calculate all distortions closer to CE as if at this point
  70 const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
  71 const Double_t AliTPCCorrection::fgkGG       =-70.0;     // Gating Grid voltage (volts)
  72
  73
  74 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
  75 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] =  {
  76 84.0,   84.5,   85.0,   85.5,   86.0,  87.0,    88.0,
  77 90.0,   92.0,   94.0,   96.0,   98.0,  100.0,  102.0,  104.0,  106.0,  108.0,
  78 110.0,  112.0,  114.0,  116.0,  118.0,  120.0,  122.0,  124.0,  126.0,  128.0,
  79 130.0,  132.0,  134.0,  136.0,  138.0,  140.0,  142.0,  144.0,  146.0,  148.0,
  80 150.0,  152.0,  154.0,  156.0,  158.0,  160.0,  162.0,  164.0,  166.0,  168.0,
  81 170.0,  172.0,  174.0,  176.0,  178.0,  180.0,  182.0,  184.0,  186.0,  188.0,
  82 190.0,  192.0,  194.0,  196.0,  198.0,  200.0,  202.0,  204.0,  206.0,  208.0,
  83 210.0,  212.0,  214.0,  216.0,  218.0,  220.0,  222.0,  224.0,  226.0,  228.0,
  84 230.0,  232.0,  234.0,  236.0,  238.0,  240.0,  242.0,  244.0,  246.0,  248.0,
  85 249.0,  249.5,  250.0,  251.5,  252.0  } ;
  86
  87 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ]     =   {
  88 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
  89 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
  90 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
  91 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
  92 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
  93 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
  94 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
  95 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
  96 -100.0,  -98.0,  -96.0,  -94.0,  -92.0,  -90.0,  -88.0,  -86.0,  -84.0,  -82.0,
  97 -80.0,  -78.0,  -76.0,  -74.0,  -72.0,  -70.0,  -68.0,  -66.0,  -64.0,  -62.0,
  98 -60.0,  -58.0,  -56.0,  -54.0,  -52.0,  -50.0,  -48.0,  -46.0,  -44.0,  -42.0,
  99 -40.0,  -38.0,  -36.0,  -34.0,  -32.0,  -30.0,  -28.0,  -26.0,  -24.0,  -22.0,
 100 -20.0,  -18.0,  -16.0,  -14.0,  -12.0,  -10.0,   -8.0,   -6.0,   -4.0,   -2.0,
 101 -1.0,   -0.5,   -0.2,   -0.1,  -0.05,   0.05,    0.1,    0.2,    0.5,    1.0,
 102  2.0,    4.0,    6.0,    8.0,   10.0,   12.0,   14.0,   16.0,   18.0,   20.0,
 103  22.0,   24.0,   26.0,   28.0,   30.0,   32.0,   34.0,   36.0,   38.0,   40.0,
 104  42.0,   44.0,   46.0,   48.0,   50.0,   52.0,   54.0,   56.0,   58.0,   60.0,
 105  62.0,   64.0,   66.0,   68.0,   70.0,   72.0,   74.0,   76.0,   78.0,   80.0,
 106  82.0,   84.0,   86.0,   88.0,   90.0,   92.0,   94.0,   96.0,   98.0,  100.0,
 107 102.0,  104.0,  106.0,  108.0,  110.0,  112.0,  114.0,  116.0,  118.0,  120.0,
 108 122.0,  124.0,  126.0,  128.0,  130.0,  132.0,  134.0,  136.0,  138.0,  140.0,
 109 142.0,  144.0,  146.0,  148.0,  150.0,  152.0,  154.0,  156.0,  158.0,  160.0,
 110 162.0,  164.0,  166.0,  168.0,  170.0,  172.0,  174.0,  176.0,  178.0,  180.0,
 111 182.0,  184.0,  186.0,  188.0,  190.0,  192.0,  194.0,  196.0,  198.0,  200.0,
 112 202.0,  204.0,  206.0,  208.0,  210.0,  212.0,  214.0,  216.0,  218.0,  220.0,
 113 222.0,  224.0,  226.0,  228.0,  230.0,  232.0,  234.0,  236.0,  238.0,  240.0,
 114 242.0,  243.0,  244.0,  245.0,  246.0,  247.0,  248.0,  248.5,  249.0,  249.5   } ;
 115
 116
 117
 118 AliTPCCorrection::AliTPCCorrection()
 119   : TNamed("correction_unity","unity"),fJLow(0),fKLow(0)
 120 {
 121   //
 122   // default constructor
 123   //
 124 }
 125
 126 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
 127   : TNamed(name,title),fJLow(0),fKLow(0)
 128 {
 129   //
 130   // default constructor, that set the name and title
 131   //
 132 }
 133
 134 AliTPCCorrection::~AliTPCCorrection() {
 135   //
 136   // virtual destructor
 137   //
 138 }
 139
 140 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
 141   //
 142   // Corrects the initial coordinates x (cartesian coordinates)
 143   // according to the given effect (inherited classes)
 144   // roc represents the TPC read out chamber (offline numbering convention)
 145   //
 146   Float_t dx[3];
 147   GetCorrection(x,roc,dx);
 148   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 149 }
 150
 151 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 152   //
 153   // Corrects the initial coordinates x (cartesian coordinates) and stores the new
 154   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 155   // roc represents the TPC read out chamber (offline numbering convention)
 156   //
 157   Float_t dx[3];
 158   GetCorrection(x,roc,dx);
 159   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 160 }
 161
 162 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
 163   //
 164   // Distorts the initial coordinates x (cartesian coordinates)
 165   // according to the given effect (inherited classes)
 166   // roc represents the TPC read out chamber (offline numbering convention)
 167   //
 168   Float_t dx[3];
 169   GetDistortion(x,roc,dx);
 170   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 171 }
 172
 173 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 174   //
 175   // Distorts the initial coordinates x (cartesian coordinates) and stores the new
 176   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 177   // roc represents the TPC read out chamber (offline numbering convention)
 178   //
 179   Float_t dx[3];
 180   GetDistortion(x,roc,dx);
 181   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 182 }
 183
 184 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
 185   //
 186   // This function delivers the correction values dx in respect to the inital coordinates x
 187   // roc represents the TPC read out chamber (offline numbering convention)
 188   // Note: The dx is overwritten by the inherited effectice class ...
 189   //
 190   for (Int_t j=0;j<3;++j) { dx[j]=0.; }
 191 }
 192
 193 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
 194   //
 195   // This function delivers the distortion values dx in respect to the inital coordinates x
 196   // roc represents the TPC read out chamber (offline numbering convention)
 197   //
 198   GetCorrection(x,roc,dx);
 199   for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
 200 }
 201
 202 void AliTPCCorrection::Init() {
 203   //
 204   // Initialization funtion (not used at the moment)
 205   //
 206 }
 207
 208 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
 209   //
 210   // Update function
 211   //
 212 }
 213
 214 void AliTPCCorrection::Print(Option_t* /*option*/) const {
 215   //
 216   // Print function to check which correction classes are used
 217   // option=="d" prints details regarding the setted magnitude
 218   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 219   //
 220   printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
 221 }
 222
 223 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t /*t1*/,Float_t /*t2*/) {
 224   //
 225   // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
 226   // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
 227   // calibration run
 228   //
 229   // SetOmegaTauT1T2(omegaTau, t1, t2);
 230 }
 231
 232 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
 233   //
 234   // Simple plot functionality.
 235   // Returns a 2d hisogram which represents the corrections in radial direction (dr)
 236   // in respect to position z within the XY plane.
 237   // The histogramm has nx times ny entries.
 238   //
 239
 240   TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
 241                      nx,-250.,250.,ny,-250.,250.);
 242   Float_t x[3],dx[3];
 243   x[2]=z;
 244   Int_t roc=z>0.?0:18; // FIXME
 245   for (Int_t iy=1;iy<=ny;++iy) {
 246     x[1]=h->GetYaxis()->GetBinCenter(iy);
 247     for (Int_t ix=1;ix<=nx;++ix) {
 248       x[0]=h->GetXaxis()->GetBinCenter(ix);
 249       GetCorrection(x,roc,dx);
 250       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 251       if (90.<=r0 && r0<=250.) {
 252         Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 253         h->SetBinContent(ix,iy,r1-r0);
 254       }
 255       else
 256         h->SetBinContent(ix,iy,0.);
 257     }
 258   }
 259   return h;
 260 }
 261
 262 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
 263   //
 264   // Simple plot functionality.
 265   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 266   // in respect to position z within the XY plane.
 267   // The histogramm has nx times ny entries.
 268   //
 269
 270   TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
 271                      nx,-250.,250.,ny,-250.,250.);
 272   Float_t x[3],dx[3];
 273   x[2]=z;
 274   Int_t roc=z>0.?0:18; // FIXME
 275   for (Int_t iy=1;iy<=ny;++iy) {
 276     x[1]=h->GetYaxis()->GetBinCenter(iy);
 277     for (Int_t ix=1;ix<=nx;++ix) {
 278       x[0]=h->GetXaxis()->GetBinCenter(ix);
 279       GetCorrection(x,roc,dx);
 280       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 281       if (90.<=r0 && r0<=250.) {
 282         Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 283         Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 284
 285         Float_t dphi=phi1-phi0;
 286         if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 287         if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 288
 289         h->SetBinContent(ix,iy,r0*dphi);
 290       }
 291       else
 292         h->SetBinContent(ix,iy,0.);
 293     }
 294   }
 295   return h;
 296 }
 297
 298 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
 299   //
 300   // Simple plot functionality.
 301   // Returns a 2d hisogram which represents the corrections in r direction (dr)
 302   // in respect to angle phi within the ZR plane.
 303   // The histogramm has nx times ny entries.
 304   //
 305   TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
 306                      nz,-250.,250.,nr,85.,250.);
 307   Float_t x[3],dx[3];
 308   for (Int_t ir=1;ir<=nr;++ir) {
 309     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 310     x[0]=radius*TMath::Cos(phi);
 311     x[1]=radius*TMath::Sin(phi);
 312     for (Int_t iz=1;iz<=nz;++iz) {
 313       x[2]=h->GetXaxis()->GetBinCenter(iz);
 314       Int_t roc=x[2]>0.?0:18; // FIXME
 315       GetCorrection(x,roc,dx);
 316       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 317       Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 318       h->SetBinContent(iz,ir,r1-r0);
 319     }
 320   }
 321   printf("SDF\n");
 322   return h;
 323
 324 }
 325
 326 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
 327   //
 328   // Simple plot functionality.
 329   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 330   // in respect to angle phi within the ZR plane.
 331   // The histogramm has nx times ny entries.
 332   //
 333   TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
 334                      nz,-250.,250.,nr,85.,250.);
 335   Float_t x[3],dx[3];
 336   for (Int_t iz=1;iz<=nz;++iz) {
 337     x[2]=h->GetXaxis()->GetBinCenter(iz);
 338     Int_t roc=x[2]>0.?0:18; // FIXME
 339     for (Int_t ir=1;ir<=nr;++ir) {
 340       Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 341       x[0]=radius*TMath::Cos(phi);
 342       x[1]=radius*TMath::Sin(phi);
 343       GetCorrection(x,roc,dx);
 344       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 345       Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 346       Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 347
 348       Float_t dphi=phi1-phi0;
 349       if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 350       if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 351
 352       h->SetBinContent(iz,ir,r0*dphi);
 353     }
 354   }
 355   return h;
 356 }
 357
 358 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
 359                                    const char *xlabel,const char *ylabel,const char *zlabel,
 360                                   Int_t nbinsx,Double_t xlow,Double_t xup,
 361                                   Int_t nbinsy,Double_t ylow,Double_t yup) {
 362   //
 363   // Helper function to create a 2d histogramm of given size
 364   //
 365
 366   TString hname=name;
 367   Int_t i=0;
 368   if (gDirectory) {
 369     while (gDirectory->FindObject(hname.Data())) {
 370       hname =name;
 371       hname+="_";
 372       hname+=i;
 373       ++i;
 374     }
 375   }
 376   TH2F *h=new TH2F(hname.Data(),title,
 377                    nbinsx,xlow,xup,
 378                    nbinsy,ylow,yup);
 379   h->GetXaxis()->SetTitle(xlabel);
 380   h->GetYaxis()->SetTitle(ylabel);
 381   h->GetZaxis()->SetTitle(zlabel);
 382   h->SetStats(0);
 383   return h;
 384 }
 385
 386
 387 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
 388
 389 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
 390                                                   const Double_t er[kNZ][kNR], Double_t &er_value )
 391 {
 392   //
 393   // Interpolate table - 2D interpolation
 394   //
 395   Double_t save_er[10] ;
 396
 397   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
 398   Search( kNR,   fgkRList,  r,   fKLow   ) ;
 399   if ( fJLow < 0 ) fJLow = 0 ;   // check if out of range
 400   if ( fKLow < 0 ) fKLow = 0 ;
 401   if ( fJLow + order  >=    kNZ - 1 ) fJLow =   kNZ - 1 - order ;
 402   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 403
 404   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
 405       save_er[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
 406   }
 407   er_value = Interpolate( &fgkZList[fJLow], save_er, order, z )   ;
 408
 409 }
 410
 411
 412 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
 413                                        const Int_t order, const Double_t x )
 414 {
 415   //
 416   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
 417   //
 418
 419   Double_t y ;
 420   if ( order == 2 ) {                // Quadratic Interpolation = 2
 421     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
 422     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
 423     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
 424   } else {                           // Linear Interpolation = 1
 425     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
 426   }
 427
 428   return (y);
 429
 430 }
 431
 432
 433 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low )
 434 {
 435   //
 436   // Search an ordered table by starting at the most recently used point
 437   //
 438
 439   Long_t middle, high ;
 440   Int_t  ascend = 0, increment = 1 ;
 441
 442   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
 443
 444   if ( low < 0 || low > n-1 ) {
 445     low = -1 ; high = n ;
 446   } else {                                            // Ordered Search phase
 447     if ( (Int_t)( x >= xArray[low] ) == ascend )  {
 448       if ( low == n-1 ) return ;
 449       high = low + 1 ;
 450       while ( (Int_t)( x >= xArray[high] ) == ascend ) {
 451         low = high ;
 452         increment *= 2 ;
 453         high = low + increment ;
 454         if ( high > n-1 )  {  high = n ; break ;  }
 455       }
 456     } else {
 457       if ( low == 0 )  {  low = -1 ;  return ;  }
 458       high = low - 1 ;
 459       while ( (Int_t)( x < xArray[low] ) == ascend ) {
 460         high = low ;
 461         increment *= 2 ;
 462         if ( increment >= high )  {  low = -1 ;  break ;  }
 463         else  low = high - increment ;
 464       }
 465     }
 466   }
 467
 468   while ( (high-low) != 1 ) {                     // Binary Search Phase
 469     middle = ( high + low ) / 2 ;
 470     if ( (Int_t)( x >= xArray[middle] ) == ascend )
 471       low = middle ;
 472     else
 473       high = middle ;
 474   }
 475
 476   if ( x == xArray[n-1] ) low = n-2 ;
 477   if ( x == xArray[0]   ) low = 0 ;
 478
 479 }
 480
 481
 482 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir,TTreeSRedirector *pcstream){
 483   //
 484   // Fit the track parameters - without and with distortion
 485   // 1. Space points in the TPC are simulated along the trajectory
 486   // 2. Space points distorted
 487   // 3. Fits  the non distorted and distroted track to the reference plane at refX
 488   // 4. For visualization and debugging  purposes the space points and tracks can be stored  in the tree - using the TTreeSRedirector functionality
 489   //
 490   // trackIn   - input track parameters
 491   // refX     - reference X to fit the track
 492   // dir      - direction - out=1 or in=-1
 493   // pcstream -  debug streamer to check the results
 494   //
 495   // see AliExternalTrackParam.h documentation:
 496   // track1.fP[0] - local y (rphi)
 497   // track1.fP[1] - z
 498   // track1.fP[2] - sinus of local inclination angle
 499   // track1.fP[3] - tangent of deep angle
 500   // track1.fP[4] - 1/pt
 501   AliTPCROC * roc = AliTPCROC::Instance();
 502   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
 503   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
 504   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 505
 506   const Double_t kMaxSnp = 0.85;
 507   const Double_t kSigmaY=0.1;
 508   const Double_t kSigmaZ=0.1;
 509   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 510
 511   AliExternalTrackParam  track(trackIn); //
 512   // generate points
 513   AliTrackPointArray pointArray0(npoints0);
 514   AliTrackPointArray pointArray1(npoints0);
 515   Double_t xyz[3];
 516   AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
 517   //
 518   // simulate the track
 519   Int_t npoints=0;
 520   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
 521   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
 522     AliTrackerBase::PropagateTrackTo(&track,radius,kMass,3,kTRUE,kMaxSnp);
 523     track.GetXYZ(xyz);
 524     AliTrackPoint pIn0;                               // space point
 525     AliTrackPoint pIn1;
 526     Int_t sector= (xyz[2]>0)? 0:18;
 527     pointArray0.GetPoint(pIn0,npoints);
 528     pointArray1.GetPoint(pIn1,npoints);
 529     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 530     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
 531     DistortPoint(distPoint, sector);
 532     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
 533     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
 534     //
 535     track.Rotate(alpha);
 536     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
 537     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
 538     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
 539     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
 540     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
 541     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
 542     pointArray0.AddPoint(npoints, &pIn0);
 543     pointArray1.AddPoint(npoints, &pIn1);
 544     npoints++;
 545     if (npoints>=npoints0) break;
 546   }
 547   //
 548   // refit track
 549   //
 550   AliExternalTrackParam *track0=0;
 551   AliExternalTrackParam *track1=0;
 552   AliTrackPoint   point1,point2,point3;
 553   if (dir==1) {  //make seed inner
 554     pointArray0.GetPoint(point1,1);
 555     pointArray0.GetPoint(point2,10);
 556     pointArray0.GetPoint(point3,20);
 557   }
 558   if (dir==-1){ //make seed outer
 559     pointArray0.GetPoint(point1,npoints-20);
 560     pointArray0.GetPoint(point2,npoints-10);
 561     pointArray0.GetPoint(point3,npoints-1);
 562   }
 563   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
 564   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
 565
 566
 567   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
 568     Int_t ipoint= (dir>0) ? ipoint: npoints-1-jpoint;
 569     //
 570     AliTrackPoint pIn0;
 571     AliTrackPoint pIn1;
 572     pointArray0.GetPoint(pIn0,ipoint);
 573     pointArray1.GetPoint(pIn1,ipoint);
 574     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
 575     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
 576     //
 577     AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,1,kFALSE,kMaxSnp);
 578     AliTrackerBase::PropagateTrackTo(track1,prot1.GetX(),kMass,1,kFALSE,kMaxSnp);
 579     track.GetXYZ(xyz);
 580     //
 581     Double_t pointPos[2]={0,0};
 582     Double_t pointCov[3]={0,0,0};
 583     pointPos[0]=prot0.GetY();//local y
 584     pointPos[1]=prot0.GetZ();//local z
 585     pointCov[0]=prot0.GetCov()[3];//simay^2
 586     pointCov[1]=prot0.GetCov()[4];//sigmayz
 587     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
 588     track0->Update(pointPos,pointCov);
 589     //
 590     pointPos[0]=prot1.GetY();//local y
 591     pointPos[1]=prot1.GetZ();//local z
 592     pointCov[0]=prot1.GetCov()[3];//simay^2
 593     pointCov[1]=prot1.GetCov()[4];//sigmayz
 594     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
 595     track1->Update(pointPos,pointCov);
 596   }
 597
 598   AliTrackerBase::PropagateTrackTo(track0,refX,kMass,2.,kTRUE,kMaxSnp);
 599   track1->Rotate(track0->GetAlpha());
 600   track1->PropagateTo(track0->GetX(),AliTrackerBase::GetBz());
 601
 602   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
 603     "point0.="<<&pointArray0<<   //  points
 604     "point1.="<<&pointArray1<<   //  distorted points
 605     "trackIn.="<<&track<<       //  original track
 606     "track0.="<<track0<<         //  fitted track
 607     "track1.="<<track1<<         //  fitted distorted track
 608     "\n";
 609   new(&trackIn) AliExternalTrackParam(*track0);
 610   delete track0;
 611   return track1;
 612 }
 613
 614
 615
 616
 617
 618 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
 619   //
 620   // create the distortion tree on a mesh with granularity given by step
 621   // return the tree with distortions at given position
 622   // Map is created on the mesh with given step size
 623   //
 624   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
 625   Float_t xyz[3];
 626   for (Double_t x= -250; x<250; x+=step){
 627     for (Double_t y= -250; y<250; y+=step){
 628       Double_t r    = TMath::Sqrt(x*x+y*y);
 629       if (r<80) continue;
 630       if (r>250) continue;
 631       for (Double_t z= -250; z<250; z+=step){
 632         Int_t roc=(z>0)?0:18;
 633         xyz[0]=x;
 634         xyz[1]=y;
 635         xyz[2]=z;
 636         Double_t phi  = TMath::ATan2(y,x);
 637         DistortPoint(xyz,roc);
 638         Double_t r1    = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 639         Double_t phi1  = TMath::ATan2(xyz[1],xyz[0]);
 640         if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
 641         if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
 642         Double_t dx = xyz[0]-x;
 643         Double_t dy = xyz[1]-y;
 644         Double_t dz = xyz[2]-z;
 645         Double_t dr=r1-r;
 646         Double_t drphi=(phi1-phi)*r;
 647         (*pcstream)<<"distortion"<<
 648           "x="<<x<<           // original position
 649           "y="<<y<<
 650           "z="<<z<<
 651           "r="<<r<<
 652           "phi="<<phi<<
 653           "x1="<<xyz[0]<<      // distorted position
 654           "y1="<<xyz[1]<<
 655           "z1="<<xyz[2]<<
 656           "r1="<<r1<<
 657           "phi1="<<phi1<<
 658           //
 659           "dx="<<dx<<          // delta position
 660           "dy="<<dy<<
 661           "dz="<<dz<<
 662           "dr="<<dr<<
 663           "drphi="<<drphi<<
 664           "\n";
 665       }
 666     }
 667   }
 668   delete pcstream;
 669   TFile f(Form("correction%s.root",GetName()));
 670   TTree * tree = (TTree*)f.Get("distortion");
 671   TTree * tree2= tree->CopyTree("1");
 672   tree2->SetName(Form("dist%s",GetName()));
 673   tree2->SetDirectory(0);
 674   delete tree;
 675   return tree2;
 676 }
 677
 678
 679
 680
 681 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, TObjArray * corrArray, Int_t step, Bool_t debug ){
 682   //
 683   // Make a fit tree:
 684   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
 685   // calculates partial distortions
 686   // Partial distortion is stored in the resulting tree
 687   // Output is storred in the file distortion_<dettype>_<partype>.root
 688   // Partial  distortion is stored with the name given by correction name
 689   //
 690   //
 691   // Parameters of function:
 692   // input     - input tree
 693   // dtype     - distortion type 0 - ITSTPC,  1 -TPCTRD, 2 - TPCvertex
 694   // ppype     - parameter type
 695   // corrArray - array with partial corrections
 696   // step      - skipe entries  - if 1 all entries processed - it is slow
 697   // debug     0 if debug on also space points dumped - it is slow
 698   const Int_t kMinEntries=50;
 699   Double_t phi,theta, snp, mean,rms, entries;
 700   tinput->SetBranchAddress("theta",&theta);
 701   tinput->SetBranchAddress("phi", &phi);
 702   tinput->SetBranchAddress("snp",&snp);
 703   tinput->SetBranchAddress("mean",&mean);
 704   tinput->SetBranchAddress("rms",&rms);
 705   tinput->SetBranchAddress("entries",&entries);
 706   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
 707   //
 708   Int_t nentries=tinput->GetEntries();
 709   Int_t ncorr=corrArray->GetEntries();
 710   Double_t corrections[100]; //
 711   Double_t tPar[5];
 712   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 713   Double_t refX=0;
 714   Int_t dir=0;
 715   if (dtype==0) {refX=85; dir=-1;}
 716   if (dtype==1) {refX=245; dir=1;}
 717   if (dtype==2) {refX=0; dir=-1;}
 718   //
 719   for (Int_t ientry=0; ientry<nentries; ientry+=step){
 720     tinput->GetEntry(ientry);
 721     tPar[0]=0;
 722     tPar[1]=theta*refX;
 723     tPar[2]=snp;
 724     tPar[3]=theta;
 725     tPar[4]=0.00001;  // should be calculated - non equal to 0
 726     cout<<endl<<endl;
 727     cout<<"Entry\n\n"<<ientry<<endl;
 728     cout<<"dtype="<<dtype<<   // detector match type
 729       "ptype="<<ptype<<   // parameter type
 730       "theta="<<theta<<   // theta
 731       "phi="<<phi<<       // phi
 732       "snp="<<phi<<       // snp
 733       "mean="<<mean<<     // mean dist value
 734       "rms="<<rms<<       // rms
 735       "entries="<<entries<<endl; // number of entries in bin
 736
 737     if (TMath::Abs(snp)>0.251) continue;
 738     (*pcstream)<<"fit"<<
 739       "dtype="<<dtype<<   // detector match type
 740       "ptype="<<ptype<<   // parameter type
 741       "theta="<<theta<<   // theta
 742       "phi="<<phi<<       // phi
 743       "snp="<<snp<<       // snp
 744       "mean="<<mean<<     // mean dist value
 745       "rms="<<rms<<       // rms
 746       "entries="<<entries;// number of entries in bin
 747     //
 748     for (Int_t icorr=0; icorr<ncorr; icorr++) {
 749       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
 750       corrections[icorr]=0;
 751       if (entries>kMinEntries){
 752         AliExternalTrackParam trackIn(refX,phi,tPar,cov);
 753         AliExternalTrackParam *trackOut = 0;
 754         if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
 755         if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
 756         corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
 757         delete trackOut;
 758       }
 759       (*pcstream)<<"fit"<<
 760         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
 761     }
 762     (*pcstream)<<"fit"<<"\n";
 763   }
 764   delete pcstream;
 765 }
 766
 767
 768
 769
 770 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
 771   //
 772   // Store object in the OCDB
 773   // By default the object is stored in the current directory
 774   // default comment consit of user name and the date
 775   //
 776   TString ocdbStorage="";
 777   ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
 778   AliCDBMetaData *metaData= new AliCDBMetaData();
 779   metaData->SetObjectClassName("AliTPCCorrection");
 780   metaData->SetResponsible("Marian Ivanov");
 781   metaData->SetBeamPeriod(1);
 782   metaData->SetAliRootVersion("05-25-01"); //root version
 783   TString userName=gSystem->GetFromPipe("echo $USER");
 784   TString date=gSystem->GetFromPipe("date");
 785
 786   if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
 787   if (comment) metaData->SetComment(comment);
 788   AliCDBId* id1=NULL;
 789   id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
 790   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
 791   gStorage->Put(this, (*id1), metaData);
 792 }
 793