TPC/AliTPCCorrection.cxx

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