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), fT1(1), fT2(1)
 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), fT1(1), fT2(1)
 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   fT1=t1;
 230   fT2=t2;
 231   //SetOmegaTauT1T2(omegaTau, t1, t2);
 232 }
 233
 234 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
 235   //
 236   // Simple plot functionality.
 237   // Returns a 2d hisogram which represents the corrections in radial direction (dr)
 238   // in respect to position z within the XY plane.
 239   // The histogramm has nx times ny entries.
 240   //
 241
 242   TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
 243                      nx,-250.,250.,ny,-250.,250.);
 244   Float_t x[3],dx[3];
 245   x[2]=z;
 246   Int_t roc=z>0.?0:18; // FIXME
 247   for (Int_t iy=1;iy<=ny;++iy) {
 248     x[1]=h->GetYaxis()->GetBinCenter(iy);
 249     for (Int_t ix=1;ix<=nx;++ix) {
 250       x[0]=h->GetXaxis()->GetBinCenter(ix);
 251       GetCorrection(x,roc,dx);
 252       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 253       if (90.<=r0 && r0<=250.) {
 254         Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 255         h->SetBinContent(ix,iy,r1-r0);
 256       }
 257       else
 258         h->SetBinContent(ix,iy,0.);
 259     }
 260   }
 261   return h;
 262 }
 263
 264 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
 265   //
 266   // Simple plot functionality.
 267   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 268   // in respect to position z within the XY plane.
 269   // The histogramm has nx times ny entries.
 270   //
 271
 272   TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
 273                      nx,-250.,250.,ny,-250.,250.);
 274   Float_t x[3],dx[3];
 275   x[2]=z;
 276   Int_t roc=z>0.?0:18; // FIXME
 277   for (Int_t iy=1;iy<=ny;++iy) {
 278     x[1]=h->GetYaxis()->GetBinCenter(iy);
 279     for (Int_t ix=1;ix<=nx;++ix) {
 280       x[0]=h->GetXaxis()->GetBinCenter(ix);
 281       GetCorrection(x,roc,dx);
 282       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 283       if (90.<=r0 && r0<=250.) {
 284         Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 285         Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 286
 287         Float_t dphi=phi1-phi0;
 288         if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 289         if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 290
 291         h->SetBinContent(ix,iy,r0*dphi);
 292       }
 293       else
 294         h->SetBinContent(ix,iy,0.);
 295     }
 296   }
 297   return h;
 298 }
 299
 300 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
 301   //
 302   // Simple plot functionality.
 303   // Returns a 2d hisogram which represents the corrections in r direction (dr)
 304   // in respect to angle phi within the ZR plane.
 305   // The histogramm has nx times ny entries.
 306   //
 307   TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
 308                      nz,-250.,250.,nr,85.,250.);
 309   Float_t x[3],dx[3];
 310   for (Int_t ir=1;ir<=nr;++ir) {
 311     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 312     x[0]=radius*TMath::Cos(phi);
 313     x[1]=radius*TMath::Sin(phi);
 314     for (Int_t iz=1;iz<=nz;++iz) {
 315       x[2]=h->GetXaxis()->GetBinCenter(iz);
 316       Int_t roc=x[2]>0.?0:18; // FIXME
 317       GetCorrection(x,roc,dx);
 318       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 319       Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 320       h->SetBinContent(iz,ir,r1-r0);
 321     }
 322   }
 323   printf("SDF\n");
 324   return h;
 325
 326 }
 327
 328 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
 329   //
 330   // Simple plot functionality.
 331   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 332   // in respect to angle phi within the ZR plane.
 333   // The histogramm has nx times ny entries.
 334   //
 335   TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
 336                      nz,-250.,250.,nr,85.,250.);
 337   Float_t x[3],dx[3];
 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     for (Int_t ir=1;ir<=nr;++ir) {
 342       Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 343       x[0]=radius*TMath::Cos(phi);
 344       x[1]=radius*TMath::Sin(phi);
 345       GetCorrection(x,roc,dx);
 346       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 347       Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 348       Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 349
 350       Float_t dphi=phi1-phi0;
 351       if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 352       if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 353
 354       h->SetBinContent(iz,ir,r0*dphi);
 355     }
 356   }
 357   return h;
 358 }
 359
 360 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
 361                                    const char *xlabel,const char *ylabel,const char *zlabel,
 362                                   Int_t nbinsx,Double_t xlow,Double_t xup,
 363                                   Int_t nbinsy,Double_t ylow,Double_t yup) {
 364   //
 365   // Helper function to create a 2d histogramm of given size
 366   //
 367
 368   TString hname=name;
 369   Int_t i=0;
 370   if (gDirectory) {
 371     while (gDirectory->FindObject(hname.Data())) {
 372       hname =name;
 373       hname+="_";
 374       hname+=i;
 375       ++i;
 376     }
 377   }
 378   TH2F *h=new TH2F(hname.Data(),title,
 379                    nbinsx,xlow,xup,
 380                    nbinsy,ylow,yup);
 381   h->GetXaxis()->SetTitle(xlabel);
 382   h->GetYaxis()->SetTitle(ylabel);
 383   h->GetZaxis()->SetTitle(zlabel);
 384   h->SetStats(0);
 385   return h;
 386 }
 387
 388
 389 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
 390
 391 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
 392                                                   const Double_t er[kNZ][kNR], Double_t &er_value )
 393 {
 394   //
 395   // Interpolate table - 2D interpolation
 396   //
 397   Double_t save_er[10] ;
 398
 399   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
 400   Search( kNR,   fgkRList,  r,   fKLow   ) ;
 401   if ( fJLow < 0 ) fJLow = 0 ;   // check if out of range
 402   if ( fKLow < 0 ) fKLow = 0 ;
 403   if ( fJLow + order  >=    kNZ - 1 ) fJLow =   kNZ - 1 - order ;
 404   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 405
 406   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
 407       save_er[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
 408   }
 409   er_value = Interpolate( &fgkZList[fJLow], save_er, order, z )   ;
 410
 411 }
 412
 413
 414 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
 415                                        const Int_t order, const Double_t x )
 416 {
 417   //
 418   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
 419   //
 420
 421   Double_t y ;
 422   if ( order == 2 ) {                // Quadratic Interpolation = 2
 423     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
 424     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
 425     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
 426   } else {                           // Linear Interpolation = 1
 427     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
 428   }
 429
 430   return (y);
 431
 432 }
 433
 434
 435 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low )
 436 {
 437   //
 438   // Search an ordered table by starting at the most recently used point
 439   //
 440
 441   Long_t middle, high ;
 442   Int_t  ascend = 0, increment = 1 ;
 443
 444   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
 445
 446   if ( low < 0 || low > n-1 ) {
 447     low = -1 ; high = n ;
 448   } else {                                            // Ordered Search phase
 449     if ( (Int_t)( x >= xArray[low] ) == ascend )  {
 450       if ( low == n-1 ) return ;
 451       high = low + 1 ;
 452       while ( (Int_t)( x >= xArray[high] ) == ascend ) {
 453         low = high ;
 454         increment *= 2 ;
 455         high = low + increment ;
 456         if ( high > n-1 )  {  high = n ; break ;  }
 457       }
 458     } else {
 459       if ( low == 0 )  {  low = -1 ;  return ;  }
 460       high = low - 1 ;
 461       while ( (Int_t)( x < xArray[low] ) == ascend ) {
 462         high = low ;
 463         increment *= 2 ;
 464         if ( increment >= high )  {  low = -1 ;  break ;  }
 465         else  low = high - increment ;
 466       }
 467     }
 468   }
 469
 470   while ( (high-low) != 1 ) {                     // Binary Search Phase
 471     middle = ( high + low ) / 2 ;
 472     if ( (Int_t)( x >= xArray[middle] ) == ascend )
 473       low = middle ;
 474     else
 475       high = middle ;
 476   }
 477
 478   if ( x == xArray[n-1] ) low = n-2 ;
 479   if ( x == xArray[0]   ) low = 0 ;
 480
 481 }
 482
 483
 484 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir,TTreeSRedirector *pcstream){
 485   //
 486   // Fit the track parameters - without and with distortion
 487   // 1. Space points in the TPC are simulated along the trajectory
 488   // 2. Space points distorted
 489   // 3. Fits  the non distorted and distroted track to the reference plane at refX
 490   // 4. For visualization and debugging  purposes the space points and tracks can be stored  in the tree - using the TTreeSRedirector functionality
 491   //
 492   // trackIn   - input track parameters
 493   // refX     - reference X to fit the track
 494   // dir      - direction - out=1 or in=-1
 495   // pcstream -  debug streamer to check the results
 496   //
 497   // see AliExternalTrackParam.h documentation:
 498   // track1.fP[0] - local y (rphi)
 499   // track1.fP[1] - z
 500   // track1.fP[2] - sinus of local inclination angle
 501   // track1.fP[3] - tangent of deep angle
 502   // track1.fP[4] - 1/pt
 503   AliTPCROC * roc = AliTPCROC::Instance();
 504   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
 505   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
 506   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 507
 508   const Double_t kMaxSnp = 0.85;
 509   const Double_t kSigmaY=0.1;
 510   const Double_t kSigmaZ=0.1;
 511   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 512
 513   AliExternalTrackParam  track(trackIn); //
 514   // generate points
 515   AliTrackPointArray pointArray0(npoints0);
 516   AliTrackPointArray pointArray1(npoints0);
 517   Double_t xyz[3];
 518   AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
 519   //
 520   // simulate the track
 521   Int_t npoints=0;
 522   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
 523   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
 524     AliTrackerBase::PropagateTrackTo(&track,radius,kMass,3,kTRUE,kMaxSnp);
 525     track.GetXYZ(xyz);
 526     AliTrackPoint pIn0;                               // space point
 527     AliTrackPoint pIn1;
 528     Int_t sector= (xyz[2]>0)? 0:18;
 529     pointArray0.GetPoint(pIn0,npoints);
 530     pointArray1.GetPoint(pIn1,npoints);
 531     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 532     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
 533     DistortPoint(distPoint, sector);
 534     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
 535     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
 536     //
 537     track.Rotate(alpha);
 538     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
 539     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
 540     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
 541     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
 542     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
 543     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
 544     pointArray0.AddPoint(npoints, &pIn0);
 545     pointArray1.AddPoint(npoints, &pIn1);
 546     npoints++;
 547     if (npoints>=npoints0) break;
 548   }
 549   //
 550   // refit track
 551   //
 552   AliExternalTrackParam *track0=0;
 553   AliExternalTrackParam *track1=0;
 554   AliTrackPoint   point1,point2,point3;
 555   if (dir==1) {  //make seed inner
 556     pointArray0.GetPoint(point1,1);
 557     pointArray0.GetPoint(point2,10);
 558     pointArray0.GetPoint(point3,20);
 559   }
 560   if (dir==-1){ //make seed outer
 561     pointArray0.GetPoint(point1,npoints-20);
 562     pointArray0.GetPoint(point2,npoints-10);
 563     pointArray0.GetPoint(point3,npoints-1);
 564   }
 565   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
 566   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
 567
 568
 569   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
 570     Int_t ipoint= (dir>0) ? ipoint: npoints-1-jpoint;
 571     //
 572     AliTrackPoint pIn0;
 573     AliTrackPoint pIn1;
 574     pointArray0.GetPoint(pIn0,ipoint);
 575     pointArray1.GetPoint(pIn1,ipoint);
 576     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
 577     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
 578     //
 579     AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,1,kFALSE,kMaxSnp);
 580     AliTrackerBase::PropagateTrackTo(track1,prot1.GetX(),kMass,1,kFALSE,kMaxSnp);
 581     track.GetXYZ(xyz);
 582     //
 583     Double_t pointPos[2]={0,0};
 584     Double_t pointCov[3]={0,0,0};
 585     pointPos[0]=prot0.GetY();//local y
 586     pointPos[1]=prot0.GetZ();//local z
 587     pointCov[0]=prot0.GetCov()[3];//simay^2
 588     pointCov[1]=prot0.GetCov()[4];//sigmayz
 589     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
 590     track0->Update(pointPos,pointCov);
 591     //
 592     pointPos[0]=prot1.GetY();//local y
 593     pointPos[1]=prot1.GetZ();//local z
 594     pointCov[0]=prot1.GetCov()[3];//simay^2
 595     pointCov[1]=prot1.GetCov()[4];//sigmayz
 596     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
 597     track1->Update(pointPos,pointCov);
 598   }
 599
 600   AliTrackerBase::PropagateTrackTo(track0,refX,kMass,2.,kTRUE,kMaxSnp);
 601   track1->Rotate(track0->GetAlpha());
 602   track1->PropagateTo(track0->GetX(),AliTrackerBase::GetBz());
 603
 604   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
 605     "point0.="<<&pointArray0<<   //  points
 606     "point1.="<<&pointArray1<<   //  distorted points
 607     "trackIn.="<<&track<<       //  original track
 608     "track0.="<<track0<<         //  fitted track
 609     "track1.="<<track1<<         //  fitted distorted track
 610     "\n";
 611   new(&trackIn) AliExternalTrackParam(*track0);
 612   delete track0;
 613   return track1;
 614 }
 615
 616
 617
 618
 619
 620 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
 621   //
 622   // create the distortion tree on a mesh with granularity given by step
 623   // return the tree with distortions at given position
 624   // Map is created on the mesh with given step size
 625   //
 626   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
 627   Float_t xyz[3];
 628   for (Double_t x= -250; x<250; x+=step){
 629     for (Double_t y= -250; y<250; y+=step){
 630       Double_t r    = TMath::Sqrt(x*x+y*y);
 631       if (r<80) continue;
 632       if (r>250) continue;
 633       for (Double_t z= -250; z<250; z+=step){
 634         Int_t roc=(z>0)?0:18;
 635         xyz[0]=x;
 636         xyz[1]=y;
 637         xyz[2]=z;
 638         Double_t phi  = TMath::ATan2(y,x);
 639         DistortPoint(xyz,roc);
 640         Double_t r1    = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 641         Double_t phi1  = TMath::ATan2(xyz[1],xyz[0]);
 642         if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
 643         if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
 644         Double_t dx = xyz[0]-x;
 645         Double_t dy = xyz[1]-y;
 646         Double_t dz = xyz[2]-z;
 647         Double_t dr=r1-r;
 648         Double_t drphi=(phi1-phi)*r;
 649         (*pcstream)<<"distortion"<<
 650           "x="<<x<<           // original position
 651           "y="<<y<<
 652           "z="<<z<<
 653           "r="<<r<<
 654           "phi="<<phi<<
 655           "x1="<<xyz[0]<<      // distorted position
 656           "y1="<<xyz[1]<<
 657           "z1="<<xyz[2]<<
 658           "r1="<<r1<<
 659           "phi1="<<phi1<<
 660           //
 661           "dx="<<dx<<          // delta position
 662           "dy="<<dy<<
 663           "dz="<<dz<<
 664           "dr="<<dr<<
 665           "drphi="<<drphi<<
 666           "\n";
 667       }
 668     }
 669   }
 670   delete pcstream;
 671   TFile f(Form("correction%s.root",GetName()));
 672   TTree * tree = (TTree*)f.Get("distortion");
 673   TTree * tree2= tree->CopyTree("1");
 674   tree2->SetName(Form("dist%s",GetName()));
 675   tree2->SetDirectory(0);
 676   delete tree;
 677   return tree2;
 678 }
 679
 680
 681
 682
 683 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, TObjArray * corrArray, Int_t step, Bool_t debug ){
 684   //
 685   // Make a fit tree:
 686   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
 687   // calculates partial distortions
 688   // Partial distortion is stored in the resulting tree
 689   // Output is storred in the file distortion_<dettype>_<partype>.root
 690   // Partial  distortion is stored with the name given by correction name
 691   //
 692   //
 693   // Parameters of function:
 694   // input     - input tree
 695   // dtype     - distortion type 0 - ITSTPC,  1 -TPCTRD, 2 - TPCvertex
 696   // ppype     - parameter type
 697   // corrArray - array with partial corrections
 698   // step      - skipe entries  - if 1 all entries processed - it is slow
 699   // debug     0 if debug on also space points dumped - it is slow
 700   const Int_t kMinEntries=50;
 701   Double_t phi,theta, snp, mean,rms, entries;
 702   tinput->SetBranchAddress("theta",&theta);
 703   tinput->SetBranchAddress("phi", &phi);
 704   tinput->SetBranchAddress("snp",&snp);
 705   tinput->SetBranchAddress("mean",&mean);
 706   tinput->SetBranchAddress("rms",&rms);
 707   tinput->SetBranchAddress("entries",&entries);
 708   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
 709   //
 710   Int_t nentries=tinput->GetEntries();
 711   Int_t ncorr=corrArray->GetEntries();
 712   Double_t corrections[100]; //
 713   Double_t tPar[5];
 714   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 715   Double_t refX=0;
 716   Int_t dir=0;
 717   if (dtype==0) {refX=85; dir=-1;}
 718   if (dtype==1) {refX=245; dir=1;}
 719   if (dtype==2) {refX=0; dir=-1;}
 720   //
 721   for (Int_t ientry=0; ientry<nentries; ientry+=step){
 722     tinput->GetEntry(ientry);
 723     tPar[0]=0;
 724     tPar[1]=theta*refX;
 725     tPar[2]=snp;
 726     tPar[3]=theta;
 727     tPar[4]=0.00001;  // should be calculated - non equal to 0
 728     cout<<endl<<endl;
 729     cout<<"Entry\n\n"<<ientry<<endl;
 730     cout<<"dtype="<<dtype<<   // detector match type
 731       "ptype="<<ptype<<   // parameter type
 732       "theta="<<theta<<   // theta
 733       "phi="<<phi<<       // phi
 734       "snp="<<phi<<       // snp
 735       "mean="<<mean<<     // mean dist value
 736       "rms="<<rms<<       // rms
 737       "entries="<<entries<<endl; // number of entries in bin
 738
 739     if (TMath::Abs(snp)>0.251) continue;
 740     (*pcstream)<<"fit"<<
 741       "dtype="<<dtype<<   // detector match type
 742       "ptype="<<ptype<<   // parameter type
 743       "theta="<<theta<<   // theta
 744       "phi="<<phi<<       // phi
 745       "snp="<<snp<<       // snp
 746       "mean="<<mean<<     // mean dist value
 747       "rms="<<rms<<       // rms
 748       "entries="<<entries;// number of entries in bin
 749     //
 750     for (Int_t icorr=0; icorr<ncorr; icorr++) {
 751       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
 752       corrections[icorr]=0;
 753       if (entries>kMinEntries){
 754         AliExternalTrackParam trackIn(refX,phi,tPar,cov);
 755         AliExternalTrackParam *trackOut = 0;
 756         if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
 757         if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
 758         corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
 759         delete trackOut;
 760       }
 761       (*pcstream)<<"fit"<<
 762         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
 763     }
 764     (*pcstream)<<"fit"<<"\n";
 765   }
 766   delete pcstream;
 767 }
 768
 769
 770
 771
 772 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
 773   //
 774   // Store object in the OCDB
 775   // By default the object is stored in the current directory
 776   // default comment consit of user name and the date
 777   //
 778   TString ocdbStorage="";
 779   ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
 780   AliCDBMetaData *metaData= new AliCDBMetaData();
 781   metaData->SetObjectClassName("AliTPCCorrection");
 782   metaData->SetResponsible("Marian Ivanov");
 783   metaData->SetBeamPeriod(1);
 784   metaData->SetAliRootVersion("05-25-01"); //root version
 785   TString userName=gSystem->GetFromPipe("echo $USER");
 786   TString date=gSystem->GetFromPipe("date");
 787
 788   if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
 789   if (comment) metaData->SetComment(comment);
 790   AliCDBId* id1=NULL;
 791   id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
 792   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
 793   gStorage->Put(this, (*id1), metaData);
 794 }
 795