TPC/AliTPCCorrection.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 ////////////////////////////////////////////////////////////////////////////////
  17 //                                                                            //
  18 // AliTPCCorrection class                                                     //
  19 //                                                                            //
  20 // This class provides a general framework to deal with space point           //
  21 // distortions. An correction class which inherits from here is for example   //
  22 // AliTPCExBBShape or AliTPCExBTwist                                          //
  23 //                                                                            //
  24 // General functions are (for example):                                       //
  25 //   CorrectPoint(x,roc) where x is the vector of inital positions in         //
  26 //   cartesian coordinates and roc represents the Read Out chamber number     //
  27 //   according to the offline naming convention. The vector x is overwritten  //
  28 //   with the corrected coordinates.                                          //
  29 //                                                                            //
  30 // An alternative usage would be CorrectPoint(x,roc,dx), which leaves the     //
  31 //   vector x untouched, put returns the distortions via the vector dx        //
  32 //                                                                            //
  33 // The class allows "effective Omega Tau" corrections to be shifted to the    //
  34 // single distortion classes.                                                 //
  35 //                                                                            //
  36 // Note: This class is normally used via the class AliTPCComposedCorrection   //
  37 //                                                                            //
  38 // date: 27/04/2010                                                           //
  39 // Authors: Magnus Mager, Stefan Rossegger, Jim Thomas                        //
  40 ////////////////////////////////////////////////////////////////////////////////
  41 #include "Riostream.h"
  42
  43 #include <TH2F.h>
  44 #include <TMath.h>
  45 #include <TROOT.h>
  46 #include <TTreeStream.h>
  47 #include <TTree.h>
  48 #include <TFile.h>
  49 #include <TTimeStamp.h>
  50 #include <AliCDBStorage.h>
  51 #include <AliCDBId.h>
  52 #include <AliCDBMetaData.h>
  53 #include  "TVectorD.h"
  54
  55 #include "TRandom.h"
  56 #include "AliTPCTransform.h"
  57 #include "AliTPCcalibDB.h"
  58 #include "AliTPCExB.h"
  59 #include "AliTPCCorrection.h"
  60 #include "AliTPCRecoParam.h"
  61
  62 #include  "AliExternalTrackParam.h"
  63 #include  "AliTrackPointArray.h"
  64 #include  "TDatabasePDG.h"
  65 #include  "AliTrackerBase.h"
  66 #include  "AliTPCROC.h"
  67 #include  "THnSparse.h"
  68 #include  "AliTPCLaserTrack.h"
  69
  70 #include "AliTPCCorrection.h"
  71
  72 ClassImp(AliTPCCorrection)
  73
  74 // FIXME: the following values should come from the database
  75 const Double_t AliTPCCorrection::fgkTPCZ0    =249.7;     // nominal gating grid position
  76 const Double_t AliTPCCorrection::fgkIFCRadius= 83.06;    // Mean Radius of the Inner Field Cage ( 82.43 min,  83.70 max) (cm)
  77 const Double_t AliTPCCorrection::fgkOFCRadius=254.5;     // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
  78 const Double_t AliTPCCorrection::fgkZOffSet  = 0.2;      // Offset from CE: calculate all distortions closer to CE as if at this point
  79 const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
  80 const Double_t AliTPCCorrection::fgkGG       =-70.0;     // Gating Grid voltage (volts)
  81
  82
  83 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
  84 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] =  {
  85 84.0,   84.5,   85.0,   85.5,   86.0,  87.0,    88.0,
  86 90.0,   92.0,   94.0,   96.0,   98.0,  100.0,  102.0,  104.0,  106.0,  108.0,
  87 110.0,  112.0,  114.0,  116.0,  118.0,  120.0,  122.0,  124.0,  126.0,  128.0,
  88 130.0,  132.0,  134.0,  136.0,  138.0,  140.0,  142.0,  144.0,  146.0,  148.0,
  89 150.0,  152.0,  154.0,  156.0,  158.0,  160.0,  162.0,  164.0,  166.0,  168.0,
  90 170.0,  172.0,  174.0,  176.0,  178.0,  180.0,  182.0,  184.0,  186.0,  188.0,
  91 190.0,  192.0,  194.0,  196.0,  198.0,  200.0,  202.0,  204.0,  206.0,  208.0,
  92 210.0,  212.0,  214.0,  216.0,  218.0,  220.0,  222.0,  224.0,  226.0,  228.0,
  93 230.0,  232.0,  234.0,  236.0,  238.0,  240.0,  242.0,  244.0,  246.0,  248.0,
  94 249.0,  249.5,  250.0,  251.5,  252.0  } ;
  95
  96 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ]     =   {
  97 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
  98 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
  99 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
 100 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
 101 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
 102 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
 103 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
 104 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
 105 -100.0,  -98.0,  -96.0,  -94.0,  -92.0,  -90.0,  -88.0,  -86.0,  -84.0,  -82.0,
 106 -80.0,  -78.0,  -76.0,  -74.0,  -72.0,  -70.0,  -68.0,  -66.0,  -64.0,  -62.0,
 107 -60.0,  -58.0,  -56.0,  -54.0,  -52.0,  -50.0,  -48.0,  -46.0,  -44.0,  -42.0,
 108 -40.0,  -38.0,  -36.0,  -34.0,  -32.0,  -30.0,  -28.0,  -26.0,  -24.0,  -22.0,
 109 -20.0,  -18.0,  -16.0,  -14.0,  -12.0,  -10.0,   -8.0,   -6.0,   -4.0,   -2.0,
 110 -1.0,   -0.5,   -0.2,   -0.1,  -0.05,   0.05,    0.1,    0.2,    0.5,    1.0,
 111  2.0,    4.0,    6.0,    8.0,   10.0,   12.0,   14.0,   16.0,   18.0,   20.0,
 112  22.0,   24.0,   26.0,   28.0,   30.0,   32.0,   34.0,   36.0,   38.0,   40.0,
 113  42.0,   44.0,   46.0,   48.0,   50.0,   52.0,   54.0,   56.0,   58.0,   60.0,
 114  62.0,   64.0,   66.0,   68.0,   70.0,   72.0,   74.0,   76.0,   78.0,   80.0,
 115  82.0,   84.0,   86.0,   88.0,   90.0,   92.0,   94.0,   96.0,   98.0,  100.0,
 116 102.0,  104.0,  106.0,  108.0,  110.0,  112.0,  114.0,  116.0,  118.0,  120.0,
 117 122.0,  124.0,  126.0,  128.0,  130.0,  132.0,  134.0,  136.0,  138.0,  140.0,
 118 142.0,  144.0,  146.0,  148.0,  150.0,  152.0,  154.0,  156.0,  158.0,  160.0,
 119 162.0,  164.0,  166.0,  168.0,  170.0,  172.0,  174.0,  176.0,  178.0,  180.0,
 120 182.0,  184.0,  186.0,  188.0,  190.0,  192.0,  194.0,  196.0,  198.0,  200.0,
 121 202.0,  204.0,  206.0,  208.0,  210.0,  212.0,  214.0,  216.0,  218.0,  220.0,
 122 222.0,  224.0,  226.0,  228.0,  230.0,  232.0,  234.0,  236.0,  238.0,  240.0,
 123 242.0,  243.0,  244.0,  245.0,  246.0,  247.0,  248.0,  248.5,  249.0,  249.5   } ;
 124
 125
 126
 127 AliTPCCorrection::AliTPCCorrection()
 128   : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
 129 {
 130   //
 131   // default constructor
 132   //
 133 }
 134
 135 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
 136 : TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
 137 {
 138   //
 139   // default constructor, that set the name and title
 140   //
 141 }
 142
 143 AliTPCCorrection::~AliTPCCorrection() {
 144   //
 145   // virtual destructor
 146   //
 147 }
 148
 149 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
 150   //
 151   // Corrects the initial coordinates x (cartesian coordinates)
 152   // according to the given effect (inherited classes)
 153   // roc represents the TPC read out chamber (offline numbering convention)
 154   //
 155   Float_t dx[3];
 156   GetCorrection(x,roc,dx);
 157   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 158 }
 159
 160 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 161   //
 162   // Corrects the initial coordinates x (cartesian coordinates) and stores the new
 163   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 164   // roc represents the TPC read out chamber (offline numbering convention)
 165   //
 166   Float_t dx[3];
 167   GetCorrection(x,roc,dx);
 168   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 169 }
 170
 171 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
 172   //
 173   // Distorts the initial coordinates x (cartesian coordinates)
 174   // according to the given effect (inherited classes)
 175   // roc represents the TPC read out chamber (offline numbering convention)
 176   //
 177   Float_t dx[3];
 178   GetDistortion(x,roc,dx);
 179   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 180 }
 181
 182 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 183   //
 184   // Distorts the initial coordinates x (cartesian coordinates) and stores the new
 185   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 186   // roc represents the TPC read out chamber (offline numbering convention)
 187   //
 188   Float_t dx[3];
 189   GetDistortion(x,roc,dx);
 190   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 191 }
 192
 193 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
 194   //
 195   // This function delivers the correction values dx in respect to the inital coordinates x
 196   // roc represents the TPC read out chamber (offline numbering convention)
 197   // Note: The dx is overwritten by the inherited effectice class ...
 198   //
 199   for (Int_t j=0;j<3;++j) { dx[j]=0.; }
 200 }
 201
 202 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
 203   //
 204   // This function delivers the distortion values dx in respect to the inital coordinates x
 205   // roc represents the TPC read out chamber (offline numbering convention)
 206   //
 207   GetCorrection(x,roc,dx);
 208   for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
 209 }
 210
 211 void AliTPCCorrection::Init() {
 212   //
 213   // Initialization funtion (not used at the moment)
 214   //
 215 }
 216
 217 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
 218   //
 219   // Update function
 220   //
 221 }
 222
 223 void AliTPCCorrection::Print(Option_t* /*option*/) const {
 224   //
 225   // Print function to check which correction classes are used
 226   // option=="d" prints details regarding the setted magnitude
 227   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 228   //
 229   printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
 230 }
 231
 232 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
 233   //
 234   // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
 235   // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
 236   // calibration run
 237   //
 238   fT1=t1;
 239   fT2=t2;
 240   //SetOmegaTauT1T2(omegaTau, t1, t2);
 241 }
 242
 243 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
 244   //
 245   // Simple plot functionality.
 246   // Returns a 2d hisogram which represents the corrections in radial direction (dr)
 247   // in respect to position z within the XY plane.
 248   // The histogramm has nx times ny entries.
 249   //
 250
 251   TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
 252                      nx,-250.,250.,ny,-250.,250.);
 253   Float_t x[3],dx[3];
 254   x[2]=z;
 255   Int_t roc=z>0.?0:18; // FIXME
 256   for (Int_t iy=1;iy<=ny;++iy) {
 257     x[1]=h->GetYaxis()->GetBinCenter(iy);
 258     for (Int_t ix=1;ix<=nx;++ix) {
 259       x[0]=h->GetXaxis()->GetBinCenter(ix);
 260       GetCorrection(x,roc,dx);
 261       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 262       if (90.<=r0 && r0<=250.) {
 263         Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 264         h->SetBinContent(ix,iy,r1-r0);
 265       }
 266       else
 267         h->SetBinContent(ix,iy,0.);
 268     }
 269   }
 270   return h;
 271 }
 272
 273 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
 274   //
 275   // Simple plot functionality.
 276   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 277   // in respect to position z within the XY plane.
 278   // The histogramm has nx times ny entries.
 279   //
 280
 281   TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
 282                      nx,-250.,250.,ny,-250.,250.);
 283   Float_t x[3],dx[3];
 284   x[2]=z;
 285   Int_t roc=z>0.?0:18; // FIXME
 286   for (Int_t iy=1;iy<=ny;++iy) {
 287     x[1]=h->GetYaxis()->GetBinCenter(iy);
 288     for (Int_t ix=1;ix<=nx;++ix) {
 289       x[0]=h->GetXaxis()->GetBinCenter(ix);
 290       GetCorrection(x,roc,dx);
 291       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 292       if (90.<=r0 && r0<=250.) {
 293         Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 294         Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 295
 296         Float_t dphi=phi1-phi0;
 297         if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 298         if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 299
 300         h->SetBinContent(ix,iy,r0*dphi);
 301       }
 302       else
 303         h->SetBinContent(ix,iy,0.);
 304     }
 305   }
 306   return h;
 307 }
 308
 309 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
 310   //
 311   // Simple plot functionality.
 312   // Returns a 2d hisogram which represents the corrections in r direction (dr)
 313   // in respect to angle phi within the ZR plane.
 314   // The histogramm has nx times ny entries.
 315   //
 316   TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
 317                      nz,-250.,250.,nr,85.,250.);
 318   Float_t x[3],dx[3];
 319   for (Int_t ir=1;ir<=nr;++ir) {
 320     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 321     x[0]=radius*TMath::Cos(phi);
 322     x[1]=radius*TMath::Sin(phi);
 323     for (Int_t iz=1;iz<=nz;++iz) {
 324       x[2]=h->GetXaxis()->GetBinCenter(iz);
 325       Int_t roc=x[2]>0.?0:18; // FIXME
 326       GetCorrection(x,roc,dx);
 327       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 328       Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 329       h->SetBinContent(iz,ir,r1-r0);
 330     }
 331   }
 332   printf("SDF\n");
 333   return h;
 334
 335 }
 336
 337 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
 338   //
 339   // Simple plot functionality.
 340   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 341   // in respect to angle phi within the ZR plane.
 342   // The histogramm has nx times ny entries.
 343   //
 344   TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
 345                      nz,-250.,250.,nr,85.,250.);
 346   Float_t x[3],dx[3];
 347   for (Int_t iz=1;iz<=nz;++iz) {
 348     x[2]=h->GetXaxis()->GetBinCenter(iz);
 349     Int_t roc=x[2]>0.?0:18; // FIXME
 350     for (Int_t ir=1;ir<=nr;++ir) {
 351       Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 352       x[0]=radius*TMath::Cos(phi);
 353       x[1]=radius*TMath::Sin(phi);
 354       GetCorrection(x,roc,dx);
 355       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 356       Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 357       Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 358
 359       Float_t dphi=phi1-phi0;
 360       if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 361       if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 362
 363       h->SetBinContent(iz,ir,r0*dphi);
 364     }
 365   }
 366   return h;
 367 }
 368
 369 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
 370                                    const char *xlabel,const char *ylabel,const char *zlabel,
 371                                   Int_t nbinsx,Double_t xlow,Double_t xup,
 372                                   Int_t nbinsy,Double_t ylow,Double_t yup) {
 373   //
 374   // Helper function to create a 2d histogramm of given size
 375   //
 376
 377   TString hname=name;
 378   Int_t i=0;
 379   if (gDirectory) {
 380     while (gDirectory->FindObject(hname.Data())) {
 381       hname =name;
 382       hname+="_";
 383       hname+=i;
 384       ++i;
 385     }
 386   }
 387   TH2F *h=new TH2F(hname.Data(),title,
 388                    nbinsx,xlow,xup,
 389                    nbinsy,ylow,yup);
 390   h->GetXaxis()->SetTitle(xlabel);
 391   h->GetYaxis()->SetTitle(ylabel);
 392   h->GetZaxis()->SetTitle(zlabel);
 393   h->SetStats(0);
 394   return h;
 395 }
 396
 397
 398 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
 399
 400 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
 401                                                   const Double_t er[kNZ][kNR], Double_t &erValue ) {
 402   //
 403   // Interpolate table - 2D interpolation
 404   //
 405   Double_t saveEr[10] ;
 406
 407   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
 408   Search( kNR,   fgkRList,  r,   fKLow   ) ;
 409   if ( fJLow < 0 ) fJLow = 0 ;   // check if out of range
 410   if ( fKLow < 0 ) fKLow = 0 ;
 411   if ( fJLow + order  >=    kNZ - 1 ) fJLow =   kNZ - 1 - order ;
 412   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 413
 414   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
 415       saveEr[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
 416   }
 417   erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z )   ;
 418
 419 }
 420
 421
 422 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
 423                                        const Int_t order, const Double_t x ) {
 424   //
 425   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
 426   //
 427
 428   Double_t y ;
 429   if ( order == 2 ) {                // Quadratic Interpolation = 2
 430     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
 431     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
 432     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
 433   } else {                           // Linear Interpolation = 1
 434     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
 435   }
 436
 437   return (y);
 438
 439 }
 440
 441
 442 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
 443   //
 444   // Search an ordered table by starting at the most recently used point
 445   //
 446
 447   Long_t middle, high ;
 448   Int_t  ascend = 0, increment = 1 ;
 449
 450   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
 451
 452   if ( low < 0 || low > n-1 ) {
 453     low = -1 ; high = n ;
 454   } else {                                            // Ordered Search phase
 455     if ( (Int_t)( x >= xArray[low] ) == ascend )  {
 456       if ( low == n-1 ) return ;
 457       high = low + 1 ;
 458       while ( (Int_t)( x >= xArray[high] ) == ascend ) {
 459         low = high ;
 460         increment *= 2 ;
 461         high = low + increment ;
 462         if ( high > n-1 )  {  high = n ; break ;  }
 463       }
 464     } else {
 465       if ( low == 0 )  {  low = -1 ;  return ;  }
 466       high = low - 1 ;
 467       while ( (Int_t)( x < xArray[low] ) == ascend ) {
 468         high = low ;
 469         increment *= 2 ;
 470         if ( increment >= high )  {  low = -1 ;  break ;  }
 471         else  low = high - increment ;
 472       }
 473     }
 474   }
 475
 476   while ( (high-low) != 1 ) {                     // Binary Search Phase
 477     middle = ( high + low ) / 2 ;
 478     if ( (Int_t)( x >= xArray[middle] ) == ascend )
 479       low = middle ;
 480     else
 481       high = middle ;
 482   }
 483
 484   if ( x == xArray[n-1] ) low = n-2 ;
 485   if ( x == xArray[0]   ) low = 0 ;
 486
 487 }
 488
 489
 490 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
 491   //
 492   // Fit the track parameters - without and with distortion
 493   // 1. Space points in the TPC are simulated along the trajectory
 494   // 2. Space points distorted
 495   // 3. Fits  the non distorted and distroted track to the reference plane at refX
 496   // 4. For visualization and debugging  purposes the space points and tracks can be stored  in the tree - using the TTreeSRedirector functionality
 497   //
 498   // trackIn   - input track parameters
 499   // refX     - reference X to fit the track
 500   // dir      - direction - out=1 or in=-1
 501   // pcstream -  debug streamer to check the results
 502   //
 503   // see AliExternalTrackParam.h documentation:
 504   // track1.fP[0] - local y (rphi)
 505   // track1.fP[1] - z
 506   // track1.fP[2] - sinus of local inclination angle
 507   // track1.fP[3] - tangent of deep angle
 508   // track1.fP[4] - 1/pt
 509   AliTPCROC * roc = AliTPCROC::Instance();
 510   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
 511   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
 512   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 513
 514   const Double_t kMaxSnp = 0.85;
 515   const Double_t kSigmaY=0.1;
 516   const Double_t kSigmaZ=0.1;
 517   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 518
 519   AliExternalTrackParam  track(trackIn); //
 520   // generate points
 521   AliTrackPointArray pointArray0(npoints0);
 522   AliTrackPointArray pointArray1(npoints0);
 523   Double_t xyz[3];
 524   AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
 525   //
 526   // simulate the track
 527   Int_t npoints=0;
 528   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
 529   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
 530     AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp);
 531     track.GetXYZ(xyz);
 532     xyz[0]+=gRandom->Gaus(0,0.005);
 533     xyz[1]+=gRandom->Gaus(0,0.005);
 534     xyz[2]+=gRandom->Gaus(0,0.005);
 535     AliTrackPoint pIn0;                               // space point
 536     AliTrackPoint pIn1;
 537     Int_t sector= (xyz[2]>0)? 0:18;
 538     pointArray0.GetPoint(pIn0,npoints);
 539     pointArray1.GetPoint(pIn1,npoints);
 540     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 541     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
 542     DistortPoint(distPoint, sector);
 543     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
 544     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
 545     //
 546     track.Rotate(alpha);
 547     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
 548     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
 549     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
 550     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
 551     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
 552     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
 553     pointArray0.AddPoint(npoints, &pIn0);
 554     pointArray1.AddPoint(npoints, &pIn1);
 555     npoints++;
 556     if (npoints>=npoints0) break;
 557   }
 558   if (npoints<npoints0/2) return 0;
 559   //
 560   // refit track
 561   //
 562   AliExternalTrackParam *track0=0;
 563   AliExternalTrackParam *track1=0;
 564   AliTrackPoint   point1,point2,point3;
 565   if (dir==1) {  //make seed inner
 566     pointArray0.GetPoint(point1,1);
 567     pointArray0.GetPoint(point2,30);
 568     pointArray0.GetPoint(point3,60);
 569   }
 570   if (dir==-1){ //make seed outer
 571     pointArray0.GetPoint(point1,npoints-60);
 572     pointArray0.GetPoint(point2,npoints-30);
 573     pointArray0.GetPoint(point3,npoints-1);
 574   }
 575   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
 576   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
 577
 578   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
 579     Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
 580     //
 581     AliTrackPoint pIn0;
 582     AliTrackPoint pIn1;
 583     pointArray0.GetPoint(pIn0,ipoint);
 584     pointArray1.GetPoint(pIn1,ipoint);
 585     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
 586     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
 587     //
 588     AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
 589     AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
 590     track.GetXYZ(xyz);  // distorted track also propagated to the same reference radius
 591     //
 592     Double_t pointPos[2]={0,0};
 593     Double_t pointCov[3]={0,0,0};
 594     pointPos[0]=prot0.GetY();//local y
 595     pointPos[1]=prot0.GetZ();//local z
 596     pointCov[0]=prot0.GetCov()[3];//simay^2
 597     pointCov[1]=prot0.GetCov()[4];//sigmayz
 598     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
 599     track0->Update(pointPos,pointCov);
 600     //
 601     Double_t deltaX=prot1.GetX()-prot0.GetX();   // delta X
 602     Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp()));  // deltaY due  delta X
 603     Double_t deltaZX=deltaX*track1->GetTgl();                           // deltaZ due  delta X
 604
 605     pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
 606     pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
 607     pointCov[0]=prot1.GetCov()[3];//simay^2
 608     pointCov[1]=prot1.GetCov()[4];//sigmayz
 609     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
 610     track1->Update(pointPos,pointCov);
 611   }
 612
 613   AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
 614   track1->Rotate(track0->GetAlpha());
 615   track1->PropagateTo(track0->GetX(),AliTrackerBase::GetBz());
 616
 617   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
 618     "point0.="<<&pointArray0<<   //  points
 619     "point1.="<<&pointArray1<<   //  distorted points
 620     "trackIn.="<<&track<<       //  original track
 621     "track0.="<<track0<<         //  fitted track
 622     "track1.="<<track1<<         //  fitted distorted track
 623     "\n";
 624   new(&trackIn) AliExternalTrackParam(*track0);
 625   delete track0;
 626   return track1;
 627 }
 628
 629
 630
 631
 632
 633 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
 634   //
 635   // create the distortion tree on a mesh with granularity given by step
 636   // return the tree with distortions at given position
 637   // Map is created on the mesh with given step size
 638   //
 639   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
 640   Float_t xyz[3];
 641   for (Double_t x= -250; x<250; x+=step){
 642     for (Double_t y= -250; y<250; y+=step){
 643       Double_t r    = TMath::Sqrt(x*x+y*y);
 644       if (r<80) continue;
 645       if (r>250) continue;
 646       for (Double_t z= -250; z<250; z+=step){
 647         Int_t roc=(z>0)?0:18;
 648         xyz[0]=x;
 649         xyz[1]=y;
 650         xyz[2]=z;
 651         Double_t phi  = TMath::ATan2(y,x);
 652         DistortPoint(xyz,roc);
 653         Double_t r1    = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 654         Double_t phi1  = TMath::ATan2(xyz[1],xyz[0]);
 655         if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
 656         if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
 657         Double_t dx = xyz[0]-x;
 658         Double_t dy = xyz[1]-y;
 659         Double_t dz = xyz[2]-z;
 660         Double_t dr=r1-r;
 661         Double_t drphi=(phi1-phi)*r;
 662         (*pcstream)<<"distortion"<<
 663           "x="<<x<<           // original position
 664           "y="<<y<<
 665           "z="<<z<<
 666           "r="<<r<<
 667           "phi="<<phi<<
 668           "x1="<<xyz[0]<<      // distorted position
 669           "y1="<<xyz[1]<<
 670           "z1="<<xyz[2]<<
 671           "r1="<<r1<<
 672           "phi1="<<phi1<<
 673           //
 674           "dx="<<dx<<          // delta position
 675           "dy="<<dy<<
 676           "dz="<<dz<<
 677           "dr="<<dr<<
 678           "drphi="<<drphi<<
 679           "\n";
 680       }
 681     }
 682   }
 683   delete pcstream;
 684   TFile f(Form("correction%s.root",GetName()));
 685   TTree * tree = (TTree*)f.Get("distortion");
 686   TTree * tree2= tree->CopyTree("1");
 687   tree2->SetName(Form("dist%s",GetName()));
 688   tree2->SetDirectory(0);
 689   delete tree;
 690   return tree2;
 691 }
 692
 693
 694
 695
 696 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
 697   //
 698   // Make a fit tree:
 699   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
 700   // calculates partial distortions
 701   // Partial distortion is stored in the resulting tree
 702   // Output is storred in the file distortion_<dettype>_<partype>.root
 703   // Partial  distortion is stored with the name given by correction name
 704   //
 705   //
 706   // Parameters of function:
 707   // input     - input tree
 708   // dtype     - distortion type 0 - ITSTPC,  1 -TPCTRD, 2 - TPCvertex
 709   // ppype     - parameter type
 710   // corrArray - array with partial corrections
 711   // step      - skipe entries  - if 1 all entries processed - it is slow
 712   // debug     0 if debug on also space points dumped - it is slow
 713   const Double_t kMaxSnp = 0.85;
 714   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 715   //  const Double_t kB2C=-0.299792458e-3;
 716   const Int_t kMinEntries=50;
 717   Double_t phi,theta, snp, mean,rms, entries;
 718   tinput->SetBranchAddress("theta",&theta);
 719   tinput->SetBranchAddress("phi", &phi);
 720   tinput->SetBranchAddress("snp",&snp);
 721   tinput->SetBranchAddress("mean",&mean);
 722   tinput->SetBranchAddress("rms",&rms);
 723   tinput->SetBranchAddress("entries",&entries);
 724   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
 725   //
 726   Int_t nentries=tinput->GetEntries();
 727   Int_t ncorr=corrArray->GetEntries();
 728   Double_t corrections[100]={0}; //
 729   Double_t tPar[5];
 730   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 731   Double_t refX=0;
 732   Int_t dir=0;
 733   if (dtype==0) {refX=85.; dir=-1;}
 734   if (dtype==1) {refX=275.; dir=1;}
 735   if (dtype==2) {refX=85.; dir=-1;}
 736   if (dtype==3) {refX=360.; dir=-1;}
 737   //
 738   for (Int_t ientry=0; ientry<nentries; ientry+=step){
 739     tinput->GetEntry(ientry);
 740     if (TMath::Abs(snp)>kMaxSnp) continue;
 741     tPar[0]=0;
 742     tPar[1]=theta*refX;
 743     tPar[2]=snp;
 744     tPar[3]=theta;
 745     tPar[4]=(gRandom->Rndm()-0.5)*0.02;  // should be calculated - non equal to 0
 746     Double_t bz=AliTrackerBase::GetBz();
 747     if (refX>10. && TMath::Abs(bz)>0.1 )  tPar[4]=snp/(refX*bz*kB2C*2);
 748     tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
 749     AliExternalTrackParam track(refX,phi,tPar,cov);
 750     Double_t xyz[3];
 751     track.GetXYZ(xyz);
 752     Int_t id=0;
 753     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
 754     if (ptype==4 &&bz<0) mean*=-1;  // interpret as curvature
 755     (*pcstream)<<"fit"<<
 756       "bz="<<bz<<         // magnetic filed used
 757       "dtype="<<dtype<<   // detector match type
 758       "ptype="<<ptype<<   // parameter type
 759       "theta="<<theta<<   // theta
 760       "phi="<<phi<<       // phi
 761       "snp="<<snp<<       // snp
 762       "mean="<<mean<<     // mean dist value
 763       "rms="<<rms<<       // rms
 764       "gx="<<xyz[0]<<         // global position at reference
 765       "gy="<<xyz[1]<<         // global position at reference
 766       "gz="<<xyz[2]<<         // global position at reference
 767       "dRrec="<<dRrec<<      // delta Radius in reconstruction
 768       "id="<<id<<             // track id
 769       "entries="<<entries;// number of entries in bin
 770     //
 771     for (Int_t icorr=0; icorr<ncorr; icorr++) {
 772       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
 773       corrections[icorr]=0;
 774       if (entries>kMinEntries){
 775         AliExternalTrackParam trackIn(refX,phi,tPar,cov);
 776         AliExternalTrackParam *trackOut = 0;
 777         if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
 778         if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
 779         if (dtype==0) {refX=85.; dir=-1;}
 780         if (dtype==1) {refX=275.; dir=1;}
 781         if (dtype==2) {refX=0; dir=-1;}
 782         if (dtype==3) {refX=360.; dir=-1;}
 783         //
 784         if (trackOut){
 785           AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
 786           trackOut->Rotate(trackIn.GetAlpha());
 787           trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
 788           //
 789           corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
 790           delete trackOut;
 791         }else{
 792           corrections[icorr]=0;
 793         }
 794         if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature
 795       }
 796       Double_t dRdummy=0;
 797       (*pcstream)<<"fit"<<
 798         Form("%s=",corr->GetName())<<corrections[icorr]<<   // dump correction value
 799         Form("dR%s=",corr->GetName())<<dRdummy;   // dump dummy correction value not needed for tracks
 800                                                   // for points it is neccessary
 801     }
 802     (*pcstream)<<"fit"<<"\n";
 803   }
 804   delete pcstream;
 805 }
 806
 807
 808
 809 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
 810   //
 811   // Make a laser fit tree for global minimization
 812   //
 813   const Double_t cutErrY=0.1;
 814   const Double_t cutErrZ=0.1;
 815   const Double_t kEpsilon=0.00000001;
 816   TVectorD *vecdY=0;
 817   TVectorD *vecdZ=0;
 818   TVectorD *veceY=0;
 819   TVectorD *veceZ=0;
 820   AliTPCLaserTrack *ltr=0;
 821   AliTPCLaserTrack::LoadTracks();
 822   tree->SetBranchAddress("dY.",&vecdY);
 823   tree->SetBranchAddress("dZ.",&vecdZ);
 824   tree->SetBranchAddress("eY.",&veceY);
 825   tree->SetBranchAddress("eZ.",&veceZ);
 826   tree->SetBranchAddress("LTr.",&ltr);
 827   Int_t entries= tree->GetEntries();
 828   TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
 829   Double_t bz=AliTrackerBase::GetBz();
 830   //
 831
 832   for (Int_t ientry=0; ientry<entries; ientry++){
 833     tree->GetEntry(ientry);
 834     if (!ltr->GetVecGX()){
 835       ltr->UpdatePoints();
 836     }
 837     TVectorD * delta= (itype==0)? vecdY:vecdZ;
 838     TVectorD * err= (itype==0)? veceY:veceZ;
 839
 840     for (Int_t irow=0; irow<159; irow++){
 841       Int_t nentries = 1000;
 842       if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
 843       if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
 844       Int_t dtype=4;
 845       Double_t phi   =(*ltr->GetVecPhi())[irow];
 846       Double_t theta =ltr->GetTgl();
 847       Double_t mean=delta->GetMatrixArray()[irow];
 848       Double_t gx=0,gy=0,gz=0;
 849       Double_t snp = (*ltr->GetVecP2())[irow];
 850       Double_t rms = 0.1+err->GetMatrixArray()[irow];
 851       gx = (*ltr->GetVecGX())[irow];
 852       gy = (*ltr->GetVecGY())[irow];
 853       gz = (*ltr->GetVecGZ())[irow];
 854       Int_t bundle= ltr->GetBundle();
 855       Double_t dRrec=0;
 856       //
 857       // get delta R used in reconstruction
 858       AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
 859       AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
 860       const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
 861       Double_t xyz0[3]={gx,gy,gz};
 862       Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
 863       //
 864       // old ExB correction
 865       //
 866       if(recoParam&&recoParam->GetUseExBCorrection()) {
 867         Double_t xyz1[3]={gx,gy,gz};
 868         calib->GetExB()->Correct(xyz0,xyz1);
 869         Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
 870         dRrec=oldR-newR;
 871       }
 872       if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
 873         Float_t xyz1[3]={gx,gy,gz};
 874         Int_t sector=(gz>0)?0:18;
 875         correction->CorrectPoint(xyz1, sector);
 876         Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
 877         dRrec=oldR-newR;
 878       }
 879
 880
 881       (*pcstream)<<"fit"<<
 882         "bz="<<bz<<         // magnetic filed used
 883         "dtype="<<dtype<<   // detector match type
 884         "ptype="<<itype<<   // parameter type
 885         "theta="<<theta<<   // theta
 886         "phi="<<phi<<       // phi
 887         "snp="<<snp<<       // snp
 888         "mean="<<mean<<     // mean dist value
 889         "rms="<<rms<<       // rms
 890         "gx="<<gx<<         // global position
 891         "gy="<<gy<<         // global position
 892         "gz="<<gz<<         // global position
 893         "dRrec="<<dRrec<<      // delta Radius in reconstruction
 894         "id="<<bundle<<     //bundle
 895         "entries="<<nentries;// number of entries in bin
 896       //
 897       //
 898       Double_t ky = TMath::Tan(TMath::ASin(snp));
 899       Int_t ncorr = corrArray->GetEntries();
 900       Double_t r0   = TMath::Sqrt(gx*gx+gy*gy);
 901       Double_t phi0 = TMath::ATan2(gy,gx);
 902       Double_t distortions[1000]={0};
 903       Double_t distortionsR[1000]={0};
 904       for (Int_t icorr=0; icorr<ncorr; icorr++) {
 905         AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
 906         Float_t distPoint[3]={gx,gy,gz};
 907         Int_t sector= (gz>0)? 0:18;
 908         if (r0>80){
 909           corr->DistortPoint(distPoint, sector);
 910         }
 911         Double_t value=distPoint[2]-gz;
 912         if (itype==0){
 913           Double_t r1   = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
 914           Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
 915           Double_t drphi= r0*(phi1-phi0);
 916           Double_t dr   = r1-r0;
 917           distortions[icorr]  = drphi-ky*dr;
 918           distortionsR[icorr] = dr;
 919         }
 920         (*pcstream)<<"fit"<<
 921           Form("%s=",corr->GetName())<<distortions[icorr]<<    // dump correction value
 922           Form("dR%s=",corr->GetName())<<distortionsR[icorr];   // dump correction R  value
 923       }
 924       (*pcstream)<<"fit"<<"\n";
 925     }
 926   }
 927   delete pcstream;
 928 }
 929
 930
 931
 932 void   AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
 933   //
 934   // make a distortion map out ou fthe residual histogram
 935   // Results are written to the debug streamer - pcstream
 936   // Parameters:
 937   //   his0       - input (4D) residual histogram
 938   //   pcstream   - file to write the tree
 939   //   run        - run number
 940   // marian.ivanov@cern.ch
 941   const Int_t kMinEntries=50;
 942   Int_t nbins1=his0->GetAxis(1)->GetNbins();
 943   Int_t first1=his0->GetAxis(1)->GetFirst();
 944   Int_t last1 =his0->GetAxis(1)->GetLast();
 945   //
 946   Double_t bz=AliTrackerBase::GetBz();
 947   Int_t idim[4]={0,1,2,3};
 948   for (Int_t ibin1=first1; ibin1<last1; ibin1++){   //axis 1 - theta
 949     //
 950     his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
 951     Double_t       x1= his0->GetAxis(1)->GetBinCenter(ibin1);
 952     THnSparse * his1 = his0->Projection(4,idim);  // projected histogram according range1
 953     Int_t nbins3     = his1->GetAxis(3)->GetNbins();
 954     Int_t first3     = his1->GetAxis(3)->GetFirst();
 955     Int_t last3      = his1->GetAxis(3)->GetLast();
 956     //
 957
 958     for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){   // axis 3 - local angle
 959       his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
 960       Double_t      x3= his1->GetAxis(3)->GetBinCenter(ibin3);
 961       if (ibin3<first3) {
 962         his1->GetAxis(3)->SetRangeUser(-1,1);
 963         x3=0;
 964       }
 965       THnSparse * his3= his1->Projection(4,idim);         //projected histogram according selection 3
 966       Int_t nbins2    = his3->GetAxis(2)->GetNbins();
 967       Int_t first2    = his3->GetAxis(2)->GetFirst();
 968       Int_t last2     = his3->GetAxis(2)->GetLast();
 969       //
 970       for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
 971         his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
 972         Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
 973         TH1 * hisDelta = his3->Projection(0);
 974         //
 975         Double_t entries = hisDelta->GetEntries();
 976         Double_t mean=0, rms=0;
 977         if (entries>kMinEntries){
 978           mean    = hisDelta->GetMean();
 979           rms = hisDelta->GetRMS();
 980         }
 981         (*pcstream)<<hname<<
 982           "run="<<run<<
 983           "bz="<<bz<<
 984           "theta="<<x1<<
 985           "phi="<<x2<<
 986           "snp="<<x3<<
 987           "entries="<<entries<<
 988           "mean="<<mean<<
 989           "rms="<<rms<<
 990           "\n";
 991         delete hisDelta;
 992         printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
 993       }
 994       delete his3;
 995     }
 996     delete his1;
 997   }
 998 }
 999
1000
1001
1002
1003
1004 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1005   //
1006   // Store object in the OCDB
1007   // By default the object is stored in the current directory
1008   // default comment consit of user name and the date
1009   //
1010   TString ocdbStorage="";
1011   ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1012   AliCDBMetaData *metaData= new AliCDBMetaData();
1013   metaData->SetObjectClassName("AliTPCCorrection");
1014   metaData->SetResponsible("Marian Ivanov");
1015   metaData->SetBeamPeriod(1);
1016   metaData->SetAliRootVersion("05-25-01"); //root version
1017   TString userName=gSystem->GetFromPipe("echo $USER");
1018   TString date=gSystem->GetFromPipe("date");
1019
1020   if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1021   if (comment) metaData->SetComment(comment);
1022   AliCDBId* id1=NULL;
1023   id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1024   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1025   gStorage->Put(this, (*id1), metaData);
1026 }
1027