TPC/Base/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 // Begin_Html
  19 //   <h2>  AliTPCCorrection class   </h2>
  20 //
  21 //   The AliTPCCorrection class provides a general framework to deal with space point distortions.
  22 //   An correction class which inherits from here is for example AliTPCExBBShape or AliTPCExBTwist. <br>
  23 //   General virtual functions are (for example) CorrectPoint(x,roc) where x is the vector of initial
  24 //   positions in cartesian coordinates and roc represents the read-out chamber number according to
  25 //   the offline numbering convention. The vector x is overwritten with the corrected coordinates. <br>
  26 //   An alternative usage would be CorrectPoint(x,roc,dx), which leaves the vector x untouched, but
  27 //   returns the distortions via the vector dx. <br>
  28 //   This class is normally used via the general class AliTPCComposedCorrection.
  29 //   <p>
  30 //   Furthermore, the class contains basic geometrical descriptions like field cage radii
  31 //   (fgkIFCRadius, fgkOFCRadius) and length (fgkTPCZ0) plus the voltages. Also, the definitions
  32 //   of size and widths of the fulcrums building the grid of the final look-up table, which is
  33 //   then interpolated, is defined in kNX and fgkXList).
  34 //   <p>
  35 //   All physics-model classes below are derived from this class in order to not duplicate code
  36 //   and to allow a uniform treatment of all physics models.
  37 //   <p>
  38 //   <h3> Poisson solver </h3>
  39 //   A numerical solver of the Poisson equation (relaxation technique) is implemented for 2-dimensional
  40 //   geometries (r,z) as well as for 3-dimensional problems (r,$\phi$,z). The corresponding function
  41 //   names are PoissonRelaxation?D. The relevant function arguments are the arrays of the boundary and
  42 //   initial conditions (ArrayofArrayV, ArrayofChargeDensities) as well as the grid granularity which
  43 //   is used during the calculation. These inputs can be chosen according to the needs of the physical
  44 //   effect which is supposed to be simulated. In the 3D version, different symmetry conditions can be set
  45 //   in order to reduce the calculation time (used in AliTPCFCVoltError3D).
  46 //   <p>
  47 //   <h3> Unified plotting functionality  </h3>
  48 //   Generic plot functions were implemented. They return a histogram pointer in the chosen plane of
  49 //   the TPC drift volume with a selectable grid granularity and the magnitude of the correction vector.
  50 //   For example, the function CreateHistoDZinXY(z,nx,ny) returns a 2-dimensional histogram which contains
  51 //   the longitudinal corrections $dz$ in the (x,y)-plane at the given z position with the granularity of
  52 //   nx and ny. The magnitude of the corrections is defined by the class from which this function is called.
  53 //   In the same manner, standard plots for the (r,$\phi$)-plane and for the other corrections like $dr$ and $rd\phi$ are available
  54 //   <p>
  55 //   Note: This class is normally used via the class AliTPCComposedCorrection
  56 // End_Html
  57 //
  58 // Begin_Macro(source)
  59 //   {
  60 //   gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
  61 //   TCanvas *c2 = new TCanvas("cAliTPCCorrection","cAliTPCCorrection",700,1050);  c2->Divide(2,3);
  62 //   AliTPCROCVoltError3D roc; // EXAMPLE PLOTS - SEE BELOW
  63 //   roc.SetOmegaTauT1T2(0,1,1); // B=0
  64 //   Float_t z0 = 1; // at +1 cm -> A side
  65 //   c2->cd(1); roc.CreateHistoDRinXY(1.,300,300)->Draw("cont4z");
  66 //   c2->cd(3);roc.CreateHistoDRPhiinXY(1.,300,300)->Draw("cont4z");
  67 //   c2->cd(5);roc.CreateHistoDZinXY(1.,300,300)->Draw("cont4z");
  68 //   Float_t phi0=0.5;
  69 //   c2->cd(2);roc.CreateHistoDRinZR(phi0)->Draw("surf2");
  70 //   c2->cd(4);roc.CreateHistoDRPhiinZR(phi0)->Draw("surf2");
  71 //   c2->cd(6);roc.CreateHistoDZinZR(phi0)->Draw("surf2");
  72 //   return c2;
  73 //   }
  74 // End_Macro
  75 //
  76 // Begin_Html
  77 //   <p>
  78 //   Date: 27/04/2010  <br>
  79 //   Authors: Magnus Mager, Stefan Rossegger, Jim Thomas
  80 // End_Html
  81 // _________________________________________________________________
  82
  83
  84 #include "Riostream.h"
  85
  86 #include <TH2F.h>
  87 #include <TMath.h>
  88 #include <TROOT.h>
  89 #include <TTreeStream.h>
  90 #include <TTree.h>
  91 #include <TFile.h>
  92 #include <TTimeStamp.h>
  93 #include <AliCDBStorage.h>
  94 #include <AliCDBId.h>
  95 #include <AliCDBMetaData.h>
  96 #include "TVectorD.h"
  97 #include "AliTPCParamSR.h"
  98
  99 #include "AliTPCCorrection.h"
 100 #include "AliLog.h"
 101
 102 #include "AliExternalTrackParam.h"
 103 #include "AliTrackPointArray.h"
 104 #include "TDatabasePDG.h"
 105 #include "AliTrackerBase.h"
 106 #include "AliTPCROC.h"
 107 #include "THnSparse.h"
 108
 109 #include "AliTPCLaserTrack.h"
 110 #include "AliESDVertex.h"
 111 #include "AliVertexerTracks.h"
 112 #include "TDatabasePDG.h"
 113 #include "TF1.h"
 114 #include "TRandom.h"
 115
 116 #include "TDatabasePDG.h"
 117
 118 #include "AliTPCTransform.h"
 119 #include "AliTPCcalibDB.h"
 120 #include "AliTPCExB.h"
 121
 122 #include "AliTPCRecoParam.h"
 123 #include "TLinearFitter.h"
 124
 125
 126 ClassImp(AliTPCCorrection)
 127
 128
 129 TObjArray *AliTPCCorrection::fgVisualCorrection=0;
 130 // instance of correction for visualization
 131
 132
 133 // FIXME: the following values should come from the database
 134 const Double_t AliTPCCorrection::fgkTPCZ0    = 249.7;     // nominal gating grid position
 135 const Double_t AliTPCCorrection::fgkIFCRadius=  83.5;     // radius which renders the "18 rod manifold" best -> compare calc. of Jim Thomas
 136 // compare gkIFCRadius=  83.05: Mean Radius of the Inner Field Cage ( 82.43 min,  83.70 max) (cm)
 137 const Double_t AliTPCCorrection::fgkOFCRadius= 254.5;     // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
 138 const Double_t AliTPCCorrection::fgkZOffSet  =   0.2;     // Offset from CE: calculate all distortions closer to CE as if at this point
 139 const Double_t AliTPCCorrection::fgkCathodeV = -100000.0; // Cathode Voltage (volts)
 140 const Double_t AliTPCCorrection::fgkGG       =     -70.0; // Gating Grid voltage (volts)
 141
 142 const Double_t  AliTPCCorrection::fgkdvdE = 0.0024; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
 143
 144 const Double_t AliTPCCorrection::fgkEM = -1.602176487e-19/9.10938215e-31; // charge/mass in [C/kg]
 145 const Double_t AliTPCCorrection::fgke0 = 8.854187817e-12;                 // vacuum permittivity [A·s/(V·m)]
 146
 147
 148 AliTPCCorrection::AliTPCCorrection()
 149   : TNamed("correction_unity","unity"),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
 150 {
 151   //
 152   // default constructor
 153   //
 154   if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
 155
 156   InitLookUpfulcrums();
 157
 158 }
 159
 160 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
 161 : TNamed(name,title),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
 162 {
 163   //
 164   // default constructor, that set the name and title
 165   //
 166   if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
 167
 168   InitLookUpfulcrums();
 169
 170 }
 171
 172 AliTPCCorrection::~AliTPCCorrection() {
 173   //
 174   // virtual destructor
 175   //
 176 }
 177
 178 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
 179   //
 180   // Corrects the initial coordinates x (cartesian coordinates)
 181   // according to the given effect (inherited classes)
 182   // roc represents the TPC read out chamber (offline numbering convention)
 183   //
 184   Float_t dx[3];
 185   GetCorrection(x,roc,dx);
 186   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 187 }
 188
 189 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 190   //
 191   // Corrects the initial coordinates x (cartesian coordinates) and stores the new
 192   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 193   // roc represents the TPC read out chamber (offline numbering convention)
 194   //
 195   Float_t dx[3];
 196   GetCorrection(x,roc,dx);
 197   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 198 }
 199
 200 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
 201   //
 202   // Distorts the initial coordinates x (cartesian coordinates)
 203   // according to the given effect (inherited classes)
 204   // roc represents the TPC read out chamber (offline numbering convention)
 205   //
 206   Float_t dx[3];
 207   GetDistortion(x,roc,dx);
 208   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 209 }
 210
 211 void AliTPCCorrection::DistortPointLocal(Float_t x[],const Short_t roc) {
 212   //
 213   // Distorts the initial coordinates x (cartesian coordinates)
 214   // according to the given effect (inherited classes)
 215   // roc represents the TPC read out chamber (offline numbering convention)
 216   //
 217   Float_t gxyz[3]={0,0,0};
 218   Double_t alpha = TMath::Pi()*(roc%18+0.5)/18;
 219   Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
 220   gxyz[0]=  ca*x[0]+sa*x[1];
 221   gxyz[1]= -sa*x[0]+ca*x[1];
 222   gxyz[2]= x[2];
 223   DistortPoint(gxyz,roc);
 224   x[0]=  ca*gxyz[0]-sa*gxyz[1];
 225   x[1]= +sa*gxyz[0]+ca*gxyz[1];
 226   x[2]= gxyz[2];
 227 }
 228 void AliTPCCorrection::CorrectPointLocal(Float_t x[],const Short_t roc) {
 229   //
 230   // Distorts the initial coordinates x (cartesian coordinates)
 231   // according to the given effect (inherited classes)
 232   // roc represents the TPC read out chamber (offline numbering convention)
 233   //
 234   Float_t gxyz[3]={0,0,0};
 235   Double_t alpha = TMath::Pi()*(roc%18+0.5)/18;
 236   Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
 237   gxyz[0]=  ca*x[0]+sa*x[1];
 238   gxyz[1]= -sa*x[0]+ca*x[1];
 239   gxyz[2]= x[2];
 240   CorrectPoint(gxyz,roc);
 241   x[0]=  ca*gxyz[0]-sa*gxyz[1];
 242   x[1]=  sa*gxyz[0]+ca*gxyz[1];
 243   x[2]=  gxyz[2];
 244 }
 245
 246 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
 247   //
 248   // Distorts the initial coordinates x (cartesian coordinates) and stores the new
 249   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
 250   // roc represents the TPC read out chamber (offline numbering convention)
 251   //
 252   Float_t dx[3];
 253   GetDistortion(x,roc,dx);
 254   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
 255 }
 256
 257 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
 258   //
 259   // This function delivers the correction values dx in respect to the inital coordinates x
 260   // roc represents the TPC read out chamber (offline numbering convention)
 261   // Note: The dx is overwritten by the inherited effectice class ...
 262   //
 263   for (Int_t j=0;j<3;++j) { dx[j]=0.; }
 264 }
 265
 266 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
 267   //
 268   // This function delivers the distortion values dx in respect to the inital coordinates x
 269   // roc represents the TPC read out chamber (offline numbering convention)
 270   //
 271   GetCorrection(x,roc,dx);
 272   for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
 273 }
 274
 275 void AliTPCCorrection::GetCorrectionDz(const Float_t x[],const Short_t roc,Float_t dx[], Float_t delta) {
 276   // author: marian.ivanov@cern.ch
 277   //
 278   // In this (virtual)function calculates the dx'/dz,  dy'/dz  and dz'/dz at given point (x,y,z)
 279   // Generic implementation. Better precision can be acchieved knowing the internal structure
 280   // of underlying trasnformation. Derived classes can reimplement it.
 281   // To calculate correction is fitted in small neighberhood:
 282   // (x+-delta,y+-delta,z+-delta) where delta is an argument
 283   //
 284   // Input parameters:
 285   //   x[]   - space point corrdinate
 286   //   roc   - readout chamber identifier (important e.g to do not miss the side of detector)
 287   //   delta - define the size of neighberhood
 288   // Output parameter:
 289   //   dx[] - array {dx'/dz,  dy'/dz ,  dz'/dz }
 290
 291   static TLinearFitter fitx(2,"pol1");
 292   static TLinearFitter fity(2,"pol1");
 293   static TLinearFitter fitz(2,"pol1");
 294   fitx.ClearPoints();
 295   fity.ClearPoints();
 296   fitz.ClearPoints();
 297   Int_t zmin=-2;
 298   Int_t zmax=0;
 299   if ((roc%36)>=18) {
 300     zmin=0;
 301     zmax=2;
 302   }
 303   //TODO: in principle one shuld check that x[2]+zdelta*delta does not get 'out of' bounds,
 304   //      so close to the CE it doesn't change the sign, since then the corrections will be wrong ...
 305   for (Int_t xdelta=-1; xdelta<=1; xdelta++)
 306     for (Int_t ydelta=-1; ydelta<=1; ydelta++){
 307 //       for (Int_t zdelta=-1; zdelta<=1; zdelta++){
 308 //   for (Int_t xdelta=-2; xdelta<=0; xdelta++)
 309 //     for (Int_t ydelta=-2; ydelta<=0; ydelta++){
 310       for (Int_t zdelta=zmin; zdelta<=zmax; zdelta++){
 311         //TODO: what happens if x[2] is on the A-Side, but x[2]+zdelta*delta
 312         //      will be on the C-Side?
 313         Float_t xyz[3]={x[0]+xdelta*delta, x[1]+ydelta*delta, x[2]+zdelta*delta};
 314         Float_t dxyz[3];
 315         GetCorrection(xyz,roc,dxyz);
 316         Double_t adelta=zdelta*delta;
 317         fitx.AddPoint(&adelta, dxyz[0]);
 318         fity.AddPoint(&adelta, dxyz[1]);
 319         fitz.AddPoint(&adelta, dxyz[2]);
 320       }
 321     }
 322   fitx.Eval();
 323   fity.Eval();
 324   fitz.Eval();
 325   dx[0] = fitx.GetParameter(1);
 326   dx[1] = fity.GetParameter(1);
 327   dx[2] = fitz.GetParameter(1);
 328 }
 329
 330 void AliTPCCorrection::GetDistortionDz(const Float_t x[],const Short_t roc,Float_t dx[], Float_t delta) {
 331   // author: marian.ivanov@cern.ch
 332   //
 333   // In this (virtual)function calculates the dx'/dz,  dy'/dz  and dz'/dz at given point (x,y,z)
 334   // Generic implementation. Better precision can be acchieved knowing the internal structure
 335   // of underlying trasnformation. Derived classes can reimplement it.
 336   // To calculate distortion is fitted in small neighberhood:
 337   // (x+-delta,y+-delta,z+-delta) where delta is an argument
 338   //
 339   // Input parameters:
 340   //   x[]   - space point corrdinate
 341   //   roc   - readout chamber identifier (important e.g to do not miss the side of detector)
 342   //   delta - define the size of neighberhood
 343   // Output parameter:
 344   //   dx[] - array {dx'/dz,  dy'/dz ,  dz'/dz }
 345
 346   static TLinearFitter fitx(2,"pol1");
 347   static TLinearFitter fity(2,"pol1");
 348   static TLinearFitter fitz(2,"pol1");
 349   fitx.ClearPoints();
 350   fity.ClearPoints();
 351   fitz.ClearPoints();
 352   //TODO: in principle one shuld check that x[2]+zdelta*delta does not get 'out of' bounds,
 353   //      so close to the CE it doesn't change the sign, since then the corrections will be wrong ...
 354   for (Int_t xdelta=-1; xdelta<=1; xdelta++)
 355     for (Int_t ydelta=-1; ydelta<=1; ydelta++){
 356       for (Int_t zdelta=-1; zdelta<=1; zdelta++){
 357         //TODO: what happens if x[2] is on the A-Side, but x[2]+zdelta*delta
 358         //      will be on the C-Side?
 359         //TODO: For the C-Side, does this have the correct sign?
 360         Float_t xyz[3]={x[0]+xdelta*delta, x[1]+ydelta*delta, x[2]+zdelta*delta};
 361         Float_t dxyz[3];
 362         GetDistortion(xyz,roc,dxyz);
 363         Double_t adelta=zdelta*delta;
 364         fitx.AddPoint(&adelta, dxyz[0]);
 365         fity.AddPoint(&adelta, dxyz[1]);
 366         fitz.AddPoint(&adelta, dxyz[2]);
 367       }
 368     }
 369     fitx.Eval();
 370     fity.Eval();
 371     fitz.Eval();
 372     dx[0] = fitx.GetParameter(1);
 373     dx[1] = fity.GetParameter(1);
 374     dx[2] = fitz.GetParameter(1);
 375 }
 376
 377 void AliTPCCorrection::GetCorrectionIntegralDz(const Float_t x[],const Short_t roc,Float_t dx[], Float_t delta){
 378   //
 379   // Integrate 3D distortion along drift lines starting from the roc plane
 380   //   to the expected z position of the point, this assumes that dz is small
 381   //   and the error propagating to z' instead of the correct z is negligible
 382   // To define the drift lines virtual function  AliTPCCorrection::GetCorrectionDz is used
 383   //
 384   // Input parameters:
 385   //   x[]   - space point corrdinate
 386   //   roc   - readout chamber identifier (important e.g to do not miss the side of detector)
 387   //   delta - define the size of neighberhood
 388   // Output parameter:
 389   //   dx[] - array { integral(dx'/dz),  integral(dy'/dz) ,  integral(dz'/dz) }
 390
 391   Float_t zroc= ((roc%36)<18) ? fgkTPCZ0:-fgkTPCZ0;
 392   Double_t zdrift = TMath::Abs(x[2]-zroc);
 393   Int_t    nsteps = Int_t(zdrift/delta)+1;
 394   //
 395   //
 396   Float_t xyz[3]={x[0],x[1],zroc};
 397   Float_t dxyz[3]={x[0],x[1],x[2]};
 398   Short_t side=(roc/18)%2;
 399   Float_t sign=1-2*side;
 400   Double_t sumdz=0;
 401   for (Int_t i=0;i<nsteps; i++){
 402     //propagate backwards, therefore opposite signs
 403     Float_t deltaZ=delta*(-sign);
 404 //     if (xyz[2]+deltaZ>fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-fgkTPCZ0);
 405 //     if (xyz[2]-deltaZ<-fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-fgkTPCZ0);
 406     // protect again integrating through the CE
 407     if (side==0){
 408       if (xyz[2]+deltaZ<0) deltaZ=-xyz[2]+1e-20;
 409     } else {
 410       if (xyz[2]+deltaZ>0) deltaZ=xyz[2]-+1e-20;
 411     }
 412     // since at larger drift (smaller z) the corrections are larger (absolute, but negative)
 413     //  the slopes will be positive.
 414     // but since we chose deltaZ opposite sign the singn of the corretion should be fine
 415
 416     GetCorrectionDz(xyz,roc,dxyz,delta);
 417     xyz[0]+=deltaZ*dxyz[0];
 418     xyz[1]+=deltaZ*dxyz[1];
 419     xyz[2]+=deltaZ;           //
 420     sumdz+=deltaZ*dxyz[2];
 421   }
 422   //
 423   dx[0]=xyz[0]-x[0];
 424   dx[1]=xyz[1]-x[1];
 425   dx[2]=      sumdz; //TODO: is sumdz correct?
 426 }
 427
 428 void AliTPCCorrection::GetDistortionIntegralDz(const Float_t x[],const Short_t roc,Float_t dx[], Float_t delta){
 429   //
 430   // Integrate 3D distortion along drift lines
 431   // To define the drift lines virtual function  AliTPCCorrection::GetCorrectionDz is used
 432   //
 433   // Input parameters:
 434   //   x[]   - space point corrdinate
 435   //   roc   - readout chamber identifier (important e.g to do not miss the side of detector)
 436   //   delta - define the size of neighberhood
 437   // Output parameter:
 438   //   dx[] - array { integral(dx'/dz),  integral(dy'/dz) ,  integral(dz'/dz) }
 439
 440   Float_t zroc= ((roc%36)<18) ? fgkTPCZ0:-fgkTPCZ0;
 441   Double_t zdrift = TMath::Abs(x[2]-zroc);
 442   Int_t    nsteps = Int_t(zdrift/delta)+1;
 443   //
 444   //
 445   Float_t xyz[3]={x[0],x[1],x[2]};
 446   Float_t dxyz[3]={x[0],x[1],x[2]};
 447   Float_t sign=((roc%36)<18) ? 1.:-1.;
 448   Double_t sumdz=0;
 449   for (Int_t i=0;i<nsteps; i++){
 450     Float_t deltaZ=delta;
 451     if (xyz[2]+deltaZ>fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-zroc);
 452     if (xyz[2]-deltaZ<-fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-zroc);
 453     // since at larger drift (smaller z) the distortions are larger
 454     //  the slopes will be negative.
 455     // and since we are moving towards the read-out plane the deltaZ for
 456     //   weighting the dK/dz should have the opposite sign
 457     deltaZ*=sign;
 458     GetDistortionDz(xyz,roc,dxyz,delta);
 459     xyz[0]+=-deltaZ*dxyz[0];
 460     xyz[1]+=-deltaZ*dxyz[1];
 461     xyz[2]+=deltaZ;           //TODO: Should this also be corrected for the dxyz[2]
 462     sumdz+=-deltaZ*dxyz[2];
 463   }
 464   //
 465   dx[0]=xyz[0]-x[0];
 466   dx[1]=xyz[1]-x[1];
 467   dx[2]=      sumdz;  //TODO: is sumdz correct?
 468
 469 }
 470
 471
 472 void AliTPCCorrection::Init() {
 473   //
 474   // Initialization funtion (not used at the moment)
 475   //
 476 }
 477
 478 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
 479   //
 480   // Update function
 481   //
 482 }
 483
 484 void AliTPCCorrection::Print(Option_t* /*option*/) const {
 485   //
 486   // Print function to check which correction classes are used
 487   // option=="d" prints details regarding the setted magnitude
 488   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 489   //
 490   printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
 491 }
 492
 493 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
 494   //
 495   // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
 496   // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
 497   // calibration run
 498   //
 499   fT1=t1;
 500   fT2=t2;
 501   //SetOmegaTauT1T2(omegaTau, t1, t2);
 502 }
 503
 504 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
 505   //
 506   // Simple plot functionality.
 507   // Returns a 2d hisogram which represents the corrections in radial direction (dr)
 508   // in respect to position z within the XY plane.
 509   // The histogramm has nx times ny entries.
 510   //
 511   AliTPCParam* tpcparam = new AliTPCParamSR;
 512
 513   TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
 514                      nx,-250.,250.,ny,-250.,250.);
 515   Float_t x[3],dx[3];
 516   x[2]=z;
 517   Int_t roc=z>0.?0:18; // FIXME
 518   for (Int_t iy=1;iy<=ny;++iy) {
 519     x[1]=h->GetYaxis()->GetBinCenter(iy);
 520     for (Int_t ix=1;ix<=nx;++ix) {
 521       x[0]=h->GetXaxis()->GetBinCenter(ix);
 522       GetCorrection(x,roc,dx);
 523       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 524       if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
 525         Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 526         h->SetBinContent(ix,iy,r1-r0);
 527       }
 528       else
 529         h->SetBinContent(ix,iy,0.);
 530     }
 531   }
 532   delete tpcparam;
 533   return h;
 534 }
 535
 536 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
 537   //
 538   // Simple plot functionality.
 539   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 540   // in respect to position z within the XY plane.
 541   // The histogramm has nx times ny entries.
 542   //
 543
 544   AliTPCParam* tpcparam = new AliTPCParamSR;
 545
 546   TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
 547                      nx,-250.,250.,ny,-250.,250.);
 548   Float_t x[3],dx[3];
 549   x[2]=z;
 550   Int_t roc=z>0.?0:18; // FIXME
 551   for (Int_t iy=1;iy<=ny;++iy) {
 552     x[1]=h->GetYaxis()->GetBinCenter(iy);
 553     for (Int_t ix=1;ix<=nx;++ix) {
 554       x[0]=h->GetXaxis()->GetBinCenter(ix);
 555       GetCorrection(x,roc,dx);
 556       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 557       if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
 558         Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 559         Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 560
 561         Float_t dphi=phi1-phi0;
 562         if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 563         if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 564
 565         h->SetBinContent(ix,iy,r0*dphi);
 566       }
 567       else
 568         h->SetBinContent(ix,iy,0.);
 569     }
 570   }
 571   delete tpcparam;
 572   return h;
 573 }
 574
 575 TH2F* AliTPCCorrection::CreateHistoDZinXY(Float_t z,Int_t nx,Int_t ny) {
 576   //
 577   // Simple plot functionality.
 578   // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
 579   // in respect to position z within the XY plane.
 580   // The histogramm has nx times ny entries.
 581   //
 582
 583   AliTPCParam* tpcparam = new AliTPCParamSR;
 584
 585   TH2F *h=CreateTH2F("dz_xy",GetTitle(),"x [cm]","y [cm]","dz [cm]",
 586                      nx,-250.,250.,ny,-250.,250.);
 587   Float_t x[3],dx[3];
 588   x[2]=z;
 589   Int_t roc=z>0.?0:18; // FIXME
 590   for (Int_t iy=1;iy<=ny;++iy) {
 591     x[1]=h->GetYaxis()->GetBinCenter(iy);
 592     for (Int_t ix=1;ix<=nx;++ix) {
 593       x[0]=h->GetXaxis()->GetBinCenter(ix);
 594       GetCorrection(x,roc,dx);
 595       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 596       if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
 597         h->SetBinContent(ix,iy,dx[2]);
 598       }
 599       else
 600         h->SetBinContent(ix,iy,0.);
 601     }
 602   }
 603   delete tpcparam;
 604   return h;
 605 }
 606
 607 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
 608   //
 609   // Simple plot functionality.
 610   // Returns a 2d hisogram which represents the corrections in r direction (dr)
 611   // in respect to angle phi within the ZR plane.
 612   // The histogramm has nx times ny entries.
 613   //
 614   TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
 615                      nz,-250.,250.,nr,85.,250.);
 616   Float_t x[3],dx[3];
 617   for (Int_t ir=1;ir<=nr;++ir) {
 618     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 619     x[0]=radius*TMath::Cos(phi);
 620     x[1]=radius*TMath::Sin(phi);
 621     for (Int_t iz=1;iz<=nz;++iz) {
 622       x[2]=h->GetXaxis()->GetBinCenter(iz);
 623       Int_t roc=x[2]>0.?0:18; // FIXME
 624       GetCorrection(x,roc,dx);
 625       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 626       Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
 627       h->SetBinContent(iz,ir,r1-r0);
 628     }
 629   }
 630   return h;
 631
 632 }
 633
 634 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
 635   //
 636   // Simple plot functionality.
 637   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
 638   // in respect to angle phi within the ZR plane.
 639   // The histogramm has nx times ny entries.
 640   //
 641   TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
 642                      nz,-250.,250.,nr,85.,250.);
 643   Float_t x[3],dx[3];
 644   for (Int_t iz=1;iz<=nz;++iz) {
 645     x[2]=h->GetXaxis()->GetBinCenter(iz);
 646     Int_t roc=x[2]>0.?0:18; // FIXME
 647     for (Int_t ir=1;ir<=nr;++ir) {
 648       Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 649       x[0]=radius*TMath::Cos(phi);
 650       x[1]=radius*TMath::Sin(phi);
 651       GetCorrection(x,roc,dx);
 652       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
 653       Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
 654       Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
 655
 656       Float_t dphi=phi1-phi0;
 657       if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
 658       if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 659
 660       h->SetBinContent(iz,ir,r0*dphi);
 661     }
 662   }
 663   return h;
 664 }
 665
 666 TH2F* AliTPCCorrection::CreateHistoDZinZR(Float_t phi,Int_t nz,Int_t nr) {
 667   //
 668   // Simple plot functionality.
 669   // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
 670   // in respect to angle phi within the ZR plane.
 671   // The histogramm has nx times ny entries.
 672   //
 673   TH2F *h=CreateTH2F("dz_zr",GetTitle(),"z [cm]","r [cm]","dz [cm]",
 674                      nz,-250.,250.,nr,85.,250.);
 675   Float_t x[3],dx[3];
 676   for (Int_t ir=1;ir<=nr;++ir) {
 677     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
 678     x[0]=radius*TMath::Cos(phi);
 679     x[1]=radius*TMath::Sin(phi);
 680     for (Int_t iz=1;iz<=nz;++iz) {
 681       x[2]=h->GetXaxis()->GetBinCenter(iz);
 682       Int_t roc=x[2]>0.?0:18; // FIXME
 683       GetCorrection(x,roc,dx);
 684       h->SetBinContent(iz,ir,dx[2]);
 685     }
 686   }
 687   return h;
 688
 689 }
 690
 691
 692 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
 693                                    const char *xlabel,const char *ylabel,const char *zlabel,
 694                                   Int_t nbinsx,Double_t xlow,Double_t xup,
 695                                   Int_t nbinsy,Double_t ylow,Double_t yup) {
 696   //
 697   // Helper function to create a 2d histogramm of given size
 698   //
 699
 700   TString hname=name;
 701   Int_t i=0;
 702   if (gDirectory) {
 703     while (gDirectory->FindObject(hname.Data())) {
 704       hname =name;
 705       hname+="_";
 706       hname+=i;
 707       ++i;
 708     }
 709   }
 710   TH2F *h=new TH2F(hname.Data(),title,
 711                    nbinsx,xlow,xup,
 712                    nbinsy,ylow,yup);
 713   h->GetXaxis()->SetTitle(xlabel);
 714   h->GetYaxis()->SetTitle(ylabel);
 715   h->GetZaxis()->SetTitle(zlabel);
 716   h->SetStats(0);
 717   return h;
 718 }
 719
 720 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
 721
 722 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
 723                                                   const Double_t er[kNZ][kNR], Double_t &erValue ) {
 724   //
 725   // Interpolate table - 2D interpolation
 726   //
 727   Double_t saveEr[5] = {0,0,0,0,0};
 728
 729   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
 730   Search( kNR,   fgkRList,  r,   fKLow   ) ;
 731   if ( fJLow < 0 ) fJLow = 0 ;   // check if out of range
 732   if ( fKLow < 0 ) fKLow = 0 ;
 733   if ( fJLow + order  >=    kNZ - 1 ) fJLow =   kNZ - 1 - order ;
 734   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 735
 736   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
 737       saveEr[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
 738   }
 739   erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z )   ;
 740
 741 }
 742
 743 void AliTPCCorrection::Interpolate3DEdistortion( const Int_t order, const Double_t r, const Float_t phi, const Double_t z,
 744                                                  const Double_t er[kNZ][kNPhi][kNR], const Double_t ephi[kNZ][kNPhi][kNR], const Double_t ez[kNZ][kNPhi][kNR],
 745                                                  Double_t &erValue, Double_t &ephiValue, Double_t &ezValue) {
 746   //
 747   // Interpolate table - 3D interpolation
 748   //
 749
 750   Double_t saveEr[5]= {0,0,0,0,0};
 751   Double_t savedEr[5]= {0,0,0,0,0} ;
 752
 753   Double_t saveEphi[5]= {0,0,0,0,0};
 754   Double_t savedEphi[5]= {0,0,0,0,0} ;
 755
 756   Double_t saveEz[5]= {0,0,0,0,0};
 757   Double_t savedEz[5]= {0,0,0,0,0} ;
 758
 759   Search( kNZ,   fgkZList,   z,   fILow   ) ;
 760   Search( kNPhi, fgkPhiList, z,   fJLow   ) ;
 761   Search( kNR,   fgkRList,   r,   fKLow   ) ;
 762
 763   if ( fILow < 0 ) fILow = 0 ;   // check if out of range
 764   if ( fJLow < 0 ) fJLow = 0 ;
 765   if ( fKLow < 0 ) fKLow = 0 ;
 766
 767   if ( fILow + order  >=    kNZ - 1 ) fILow =   kNZ - 1 - order ;
 768   if ( fJLow + order  >=  kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
 769   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
 770
 771   for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
 772     for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
 773       saveEr[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r )   ;
 774       saveEphi[j-fJLow]   = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
 775       saveEz[j-fJLow]     = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r )   ;
 776     }
 777     savedEr[i-fILow]     = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi )   ;
 778     savedEphi[i-fILow]   = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
 779     savedEz[i-fILow]     = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi )   ;
 780   }
 781   erValue     = Interpolate( &fgkZList[fILow], savedEr, order, z )    ;
 782   ephiValue   = Interpolate( &fgkZList[fILow], savedEphi, order, z )  ;
 783   ezValue     = Interpolate( &fgkZList[fILow], savedEz, order, z )    ;
 784
 785 }
 786
 787 Double_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
 788                                               const Int_t nx,  const Int_t ny, const Double_t xv[], const Double_t yv[],
 789                                               const TMatrixD &array ) {
 790   //
 791   // Interpolate table (TMatrix format) - 2D interpolation
 792   //
 793
 794   static  Int_t jlow = 0, klow = 0 ;
 795   Double_t saveArray[5] = {0,0,0,0,0} ;
 796
 797   Search( nx,  xv,  x,   jlow  ) ;
 798   Search( ny,  yv,  y,   klow  ) ;
 799   if ( jlow < 0 ) jlow = 0 ;   // check if out of range
 800   if ( klow < 0 ) klow = 0 ;
 801   if ( jlow + order  >=    nx - 1 ) jlow =   nx - 1 - order ;
 802   if ( klow + order  >=    ny - 1 ) klow =   ny - 1 - order ;
 803
 804   for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
 805     {
 806       Double_t *ajkl = &((TMatrixD&)array)(j,klow);
 807       saveArray[j-jlow]  = Interpolate( &yv[klow], ajkl , order, y )   ;
 808     }
 809
 810   return( Interpolate( &xv[jlow], saveArray, order, x ) )   ;
 811
 812 }
 813
 814 Double_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x,   const Double_t y,   const Double_t z,
 815                                               const Int_t  nx,    const Int_t  ny,    const Int_t  nz,
 816                                               const Double_t xv[], const Double_t yv[], const Double_t zv[],
 817                                               TMatrixD **arrayofArrays ) {
 818   //
 819   // Interpolate table (TMatrix format) - 3D interpolation
 820   //
 821
 822   static  Int_t ilow = 0, jlow = 0, klow = 0 ;
 823   Double_t saveArray[5]= {0,0,0,0,0};
 824   Double_t savedArray[5]= {0,0,0,0,0} ;
 825
 826   Search( nx, xv, x, ilow   ) ;
 827   Search( ny, yv, y, jlow   ) ;
 828   Search( nz, zv, z, klow   ) ;
 829
 830   if ( ilow < 0 ) ilow = 0 ;   // check if out of range
 831   if ( jlow < 0 ) jlow = 0 ;
 832   if ( klow < 0 ) klow = 0 ;
 833
 834   if ( ilow + order  >=    nx - 1 ) ilow =   nx - 1 - order ;
 835   if ( jlow + order  >=    ny - 1 ) jlow =   ny - 1 - order ;
 836   if ( klow + order  >=    nz - 1 ) klow =   nz - 1 - order ;
 837
 838   for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
 839     {
 840       TMatrixD &table = *arrayofArrays[k] ;
 841       for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
 842         {
 843           saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y )   ;
 844         }
 845       savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x )   ;
 846     }
 847   return( Interpolate( &zv[klow], savedArray, order, z ) )   ;
 848
 849 }
 850
 851 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
 852                                        const Int_t order, const Double_t x ) {
 853   //
 854   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
 855   //
 856
 857   Double_t y ;
 858   if ( order == 2 ) {                // Quadratic Interpolation = 2
 859     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
 860     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
 861     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
 862   } else {                           // Linear Interpolation = 1
 863     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
 864   }
 865
 866   return (y);
 867
 868 }
 869
 870 Float_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
 871                                               const Int_t nx,  const Int_t ny, const Double_t xv[], const Double_t yv[],
 872                                               const TMatrixF &array ) {
 873   //
 874   // Interpolate table (TMatrix format) - 2D interpolation
 875   // Float version (in order to decrease the OCDB size)
 876   //
 877
 878   static  Int_t jlow = 0, klow = 0 ;
 879   Float_t saveArray[5] = {0.,0.,0.,0.,0.} ;
 880
 881   Search( nx,  xv,  x,   jlow  ) ;
 882   Search( ny,  yv,  y,   klow  ) ;
 883   if ( jlow < 0 ) jlow = 0 ;   // check if out of range
 884   if ( klow < 0 ) klow = 0 ;
 885   if ( jlow + order  >=    nx - 1 ) jlow =   nx - 1 - order ;
 886   if ( klow + order  >=    ny - 1 ) klow =   ny - 1 - order ;
 887
 888   for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
 889     {
 890       Float_t *ajkl = &((TMatrixF&)array)(j,klow);
 891       saveArray[j-jlow]  = Interpolate( &yv[klow], ajkl , order, y )   ;
 892     }
 893
 894   return( Interpolate( &xv[jlow], saveArray, order, x ) )   ;
 895
 896 }
 897
 898 Float_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x,   const Double_t y,   const Double_t z,
 899                                               const Int_t  nx,    const Int_t  ny,    const Int_t  nz,
 900                                               const Double_t xv[], const Double_t yv[], const Double_t zv[],
 901                                               TMatrixF **arrayofArrays ) {
 902   //
 903   // Interpolate table (TMatrix format) - 3D interpolation
 904   // Float version (in order to decrease the OCDB size)
 905   //
 906
 907   static  Int_t ilow = 0, jlow = 0, klow = 0 ;
 908   Float_t saveArray[5]= {0.,0.,0.,0.,0.};
 909   Float_t savedArray[5]= {0.,0.,0.,0.,0.} ;
 910
 911   Search( nx, xv, x, ilow   ) ;
 912   Search( ny, yv, y, jlow   ) ;
 913   Search( nz, zv, z, klow   ) ;
 914
 915   if ( ilow < 0 ) ilow = 0 ;   // check if out of range
 916   if ( jlow < 0 ) jlow = 0 ;
 917   if ( klow < 0 ) klow = 0 ;
 918
 919   if ( ilow + order  >=    nx - 1 ) ilow =   nx - 1 - order ;
 920   if ( jlow + order  >=    ny - 1 ) jlow =   ny - 1 - order ;
 921   if ( klow + order  >=    nz - 1 ) klow =   nz - 1 - order ;
 922
 923   for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
 924     {
 925       TMatrixF &table = *arrayofArrays[k] ;
 926       for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
 927         {
 928           saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y )   ;
 929         }
 930       savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x )   ;
 931     }
 932   return( Interpolate( &zv[klow], savedArray, order, z ) )   ;
 933
 934 }
 935 Float_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Float_t yArray[],
 936                                        const Int_t order, const Double_t x ) {
 937   //
 938   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
 939   // Float version (in order to decrease the OCDB size)
 940   //
 941
 942   Float_t y ;
 943   if ( order == 2 ) {                // Quadratic Interpolation = 2
 944     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
 945     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
 946     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
 947   } else {                           // Linear Interpolation = 1
 948     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
 949   }
 950
 951   return (y);
 952
 953 }
 954
 955
 956
 957 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
 958   //
 959   // Search an ordered table by starting at the most recently used point
 960   //
 961
 962   Long_t middle, high ;
 963   Int_t  ascend = 0, increment = 1 ;
 964
 965   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
 966
 967   if ( low < 0 || low > n-1 ) {
 968     low = -1 ; high = n ;
 969   } else {                                            // Ordered Search phase
 970     if ( (Int_t)( x >= xArray[low] ) == ascend )  {
 971       if ( low == n-1 ) return ;
 972       high = low + 1 ;
 973       while ( (Int_t)( x >= xArray[high] ) == ascend ) {
 974         low = high ;
 975         increment *= 2 ;
 976         high = low + increment ;
 977         if ( high > n-1 )  {  high = n ; break ;  }
 978       }
 979     } else {
 980       if ( low == 0 )  {  low = -1 ;  return ;  }
 981       high = low - 1 ;
 982       while ( (Int_t)( x < xArray[low] ) == ascend ) {
 983         high = low ;
 984         increment *= 2 ;
 985         if ( increment >= high )  {  low = -1 ;  break ;  }
 986         else  low = high - increment ;
 987       }
 988     }
 989   }
 990
 991   while ( (high-low) != 1 ) {                     // Binary Search Phase
 992     middle = ( high + low ) / 2 ;
 993     if ( (Int_t)( x >= xArray[middle] ) == ascend )
 994       low = middle ;
 995     else
 996       high = middle ;
 997   }
 998
 999   if ( x == xArray[n-1] ) low = n-2 ;
1000   if ( x == xArray[0]   ) low = 0 ;
1001
1002 }
1003
1004 void AliTPCCorrection::InitLookUpfulcrums() {
1005   //
1006   // Initialization of interpolation points - for main look up table
1007   //   (course grid in the middle, fine grid on the borders)
1008   //
1009
1010   AliTPCROC * roc = AliTPCROC::Instance();
1011   const Double_t rLow =  TMath::Floor(roc->GetPadRowRadii(0,0))-1; // first padRow plus some margin
1012
1013   // fulcrums in R
1014   fgkRList[0] = rLow;
1015   for (Int_t i = 1; i<kNR; i++) {
1016     fgkRList[i] = fgkRList[i-1] + 3.5;     // 3.5 cm spacing
1017     if (fgkRList[i]<90 ||fgkRList[i]>245)
1018        fgkRList[i] = fgkRList[i-1] + 0.5; // 0.5 cm spacing
1019     else if (fgkRList[i]<100 || fgkRList[i]>235)
1020        fgkRList[i] = fgkRList[i-1] + 1.5;  // 1.5 cm spacing
1021     else if (fgkRList[i]<120 || fgkRList[i]>225)
1022        fgkRList[i] = fgkRList[i-1] + 2.5;  // 2.5 cm spacing
1023   }
1024
1025   // fulcrums in Z
1026   fgkZList[0] = -249.5;
1027   fgkZList[kNZ-1] = 249.5;
1028   for (Int_t j = 1; j<kNZ/2; j++) {
1029     fgkZList[j] = fgkZList[j-1];
1030     if      (TMath::Abs(fgkZList[j])< 0.15)
1031       fgkZList[j] = fgkZList[j-1] + 0.09; // 0.09 cm spacing
1032     else if(TMath::Abs(fgkZList[j])< 0.6)
1033       fgkZList[j] = fgkZList[j-1] + 0.4; // 0.4 cm spacing
1034     else if      (TMath::Abs(fgkZList[j])< 2.5 || TMath::Abs(fgkZList[j])>248)
1035       fgkZList[j] = fgkZList[j-1] + 0.5; // 0.5 cm spacing
1036     else if (TMath::Abs(fgkZList[j])<10 || TMath::Abs(fgkZList[j])>235)
1037       fgkZList[j] = fgkZList[j-1] + 1.5;  // 1.5 cm spacing
1038     else if (TMath::Abs(fgkZList[j])<25 || TMath::Abs(fgkZList[j])>225)
1039       fgkZList[j] = fgkZList[j-1] + 2.5;  // 2.5 cm spacing
1040     else
1041       fgkZList[j] = fgkZList[j-1] + 4;  // 4 cm spacing
1042
1043     fgkZList[kNZ-j-1] = -fgkZList[j];
1044   }
1045
1046   // fulcrums in phi
1047   for (Int_t k = 0; k<kNPhi; k++)
1048     fgkPhiList[k] = TMath::TwoPi()*k/(kNPhi-1);
1049
1050
1051 }
1052
1053
1054 void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
1055                                            TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
1056                                            const Int_t rows, const Int_t columns, const Int_t iterations,
1057                                            const Bool_t rocDisplacement ) {
1058   //
1059   // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
1060   //
1061   // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
1062   // boundary conditions on the first and last rows, and the first and last columns.  The remainder of the
1063   // array can be blank or contain a preliminary guess at the solution.  The Charge density matrix contains
1064   // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
1065   // you wish to solve Laplaces equation however it should not contain random numbers or you will get
1066   // random numbers back as a solution.
1067   // Poisson's equation is solved by iteratively relaxing the matrix to the final solution.  In order to
1068   // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
1069   // space.  First it solves the problem on a very sparse grid by skipping rows and columns in the original
1070   // matrix.  Then it doubles the number of points and solves the problem again.  Then it doubles the
1071   // number of points and solves the problem again.  This happens several times until the maximum number
1072   // of points has been included in the array.
1073   //
1074   // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
1075   // So rows == 2**M + 1 and columns == 2**N + 1.  The number of rows and columns can be different.
1076   //
1077   // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
1078   //
1079   // Original code by Jim Thomas (STAR TPC Collaboration)
1080   //
1081
1082   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
1083
1084   const Float_t  gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
1085   const Float_t  gridSizeZ   =  fgkTPCZ0 / (columns-1) ;
1086   const Float_t  ratio       =  gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
1087
1088   TMatrixD  arrayEr(rows,columns) ;
1089   TMatrixD  arrayEz(rows,columns) ;
1090
1091   //Check that number of rows and columns is suitable for a binary expansion
1092
1093   if ( !IsPowerOfTwo(rows-1) ) {
1094     AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
1095     return;
1096   }
1097   if ( !IsPowerOfTwo(columns-1) ) {
1098     AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
1099     return;
1100   }
1101
1102   // Solve Poisson's equation in cylindrical coordinates by relaxation technique
1103   // Allow for different size grid spacing in R and Z directions
1104   // Use a binary expansion of the size of the matrix to speed up the solution of the problem
1105
1106   Int_t iOne = (rows-1)/4 ;
1107   Int_t jOne = (columns-1)/4 ;
1108   // Solve for N in 2**N, add one.
1109   Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
1110
1111   for ( Int_t count = 0 ; count < loops ; count++ ) {
1112     // Loop while the matrix expands & the resolution increases.
1113
1114     Float_t tempGridSizeR = gridSizeR * iOne ;
1115     Float_t tempRatio     = ratio * iOne * iOne / ( jOne * jOne ) ;
1116     Float_t tempFourth    = 1.0 / (2.0 + 2.0*tempRatio) ;
1117
1118     // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1119     std::vector<float> coef1(rows) ;
1120     std::vector<float> coef2(rows) ;
1121
1122     for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
1123        Float_t radius = fgkIFCRadius + i*gridSizeR ;
1124       coef1[i] = 1.0 + tempGridSizeR/(2*radius);
1125       coef2[i] = 1.0 - tempGridSizeR/(2*radius);
1126     }
1127
1128     TMatrixD sumChargeDensity(rows,columns) ;
1129
1130     for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
1131       Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
1132       for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
1133         if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1134         else {
1135           // Add up all enclosed charge density contributions within 1/2 unit in all directions
1136           Float_t weight = 0.0 ;
1137           Float_t sum    = 0.0 ;
1138           sumChargeDensity(i,j) = 0.0 ;
1139           for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
1140             for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
1141               if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
1142               else
1143                 weight = 1.0 ;
1144               // Note that this is cylindrical geometry
1145               sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1146               sum += weight*radius ;
1147             }
1148           }
1149           sumChargeDensity(i,j) /= sum ;
1150         }
1151         sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
1152        }
1153     }
1154
1155     for ( Int_t k = 1 ; k <= iterations; k++ ) {
1156       // Solve Poisson's Equation
1157       // Over-relaxation index, must be >= 1 but < 2.  Arrange for it to evolve from 2 => 1
1158       // as interations increase.
1159       Float_t overRelax   = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1160       Float_t overRelaxM1 = overRelax - 1.0 ;
1161       Float_t overRelaxtempFourth, overRelaxcoef5 ;
1162       overRelaxtempFourth = overRelax * tempFourth ;
1163       overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
1164
1165       for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
1166         for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
1167
1168           arrayV(i,j) = (   coef2[i]       *   arrayV(i-iOne,j)
1169                           + tempRatio      * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
1170                           - overRelaxcoef5 *   arrayV(i,j)
1171                           + coef1[i]       *   arrayV(i+iOne,j)
1172                           + sumChargeDensity(i,j)
1173                         ) * overRelaxtempFourth;
1174         }
1175       }
1176
1177       if ( k == iterations ) {
1178         // After full solution is achieved, copy low resolution solution into higher res array
1179         for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
1180           for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
1181
1182             if ( iOne > 1 ) {
1183               arrayV(i+iOne/2,j)                    =  ( arrayV(i+iOne,j) + arrayV(i,j)     ) / 2 ;
1184               if ( i == iOne )  arrayV(i-iOne/2,j) =  ( arrayV(0,j)       + arrayV(iOne,j) ) / 2 ;
1185             }
1186             if ( jOne > 1 ) {
1187               arrayV(i,j+jOne/2)                    =  ( arrayV(i,j+jOne) + arrayV(i,j) )     / 2 ;
1188               if ( j == jOne )  arrayV(i,j-jOne/2) =  ( arrayV(i,0)       + arrayV(i,jOne) ) / 2 ;
1189             }
1190             if ( iOne > 1 && jOne > 1 ) {
1191               arrayV(i+iOne/2,j+jOne/2) =  ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
1192               if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
1193               if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) =   ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
1194               // Note that this leaves a point at the upper left and lower right corners uninitialized.
1195               // -> Not a big deal.
1196             }
1197
1198           }
1199         }
1200       }
1201
1202     }
1203
1204     iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
1205     jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
1206
1207     sumChargeDensity.Clear();
1208   }
1209
1210   // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
1211   for ( Int_t j = 0 ; j < columns ; j++ ) {
1212     for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1213     arrayEr(0,j)      =  -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1214     arrayEr(rows-1,j) =  -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1215   }
1216
1217   // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1218   for ( Int_t i = 0 ; i < rows ; i++) {
1219     for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1220     arrayEz(i,0)         =  -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1221     arrayEz(i,columns-1) =  -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1222   }
1223
1224   for ( Int_t i = 0 ; i < rows ; i++) {
1225     // Note: go back and compare to old version of this code.  See notes below.
1226     // JT Test ... attempt to divide by real Ez not Ez to first order
1227     for ( Int_t j = 0 ; j < columns ; j++ ) {
1228       arrayEz(i,j) += ezField;
1229       // This adds back the overall Z gradient of the field (main E field component)
1230     }
1231     // Warning: (-=) assumes you are using an error potetial without the overall Field included
1232   }
1233
1234   // Integrate Er/Ez from Z to zero
1235   for ( Int_t j = 0 ; j < columns ; j++ )  {
1236     for ( Int_t i = 0 ; i < rows ; i++ ) {
1237
1238       Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
1239       arrayErOverEz(i,j) = 0.0 ;
1240       arrayDeltaEz(i,j) = 0.0 ;
1241
1242       for ( Int_t k = j ; k < columns ; k++ ) {
1243         arrayErOverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
1244         arrayDeltaEz(i,j)   +=  index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
1245         if ( index != 4 )  index = 4; else index = 2 ;
1246       }
1247       if ( index == 4 ) {
1248         arrayErOverEz(i,j)  -=  (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
1249         arrayDeltaEz(i,j)   -=  (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
1250       }
1251       if ( index == 2 ) {
1252         arrayErOverEz(i,j)  +=  (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
1253                                                     -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
1254         arrayDeltaEz(i,j)   +=  (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
1255                                                     -2.5*(arrayEz(i,columns-1)-ezField));
1256       }
1257       if ( j == columns-2 ) {
1258         arrayErOverEz(i,j) =  (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
1259                                                   +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
1260         arrayDeltaEz(i,j)  =  (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
1261                                                   +1.5*(arrayEz(i,columns-1)-ezField) ) ;
1262       }
1263       if ( j == columns-1 ) {
1264         arrayErOverEz(i,j) =  0.0 ;
1265         arrayDeltaEz(i,j)  =  0.0 ;
1266       }
1267     }
1268   }
1269
1270   // calculate z distortion from the integrated Delta Ez residuals
1271   // and include the aquivalence (Volt to cm) of the ROC shift !!
1272
1273   for ( Int_t j = 0 ; j < columns ; j++ )  {
1274     for ( Int_t i = 0 ; i < rows ; i++ ) {
1275
1276       // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1277       arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
1278
1279       // ROC Potential in cm aquivalent
1280       Double_t dzROCShift =  arrayV(i, columns -1)/ezField;
1281       if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift;  // add the ROC misaligment
1282
1283     }
1284   }
1285
1286   arrayEr.Clear();
1287   arrayEz.Clear();
1288
1289 }
1290
1291 void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
1292                     TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
1293                     const Int_t rows, const Int_t columns,  const Int_t phislices,
1294                     const Float_t deltaphi, const Int_t iterations, const Int_t symmetry,
1295                     Bool_t rocDisplacement  ) {
1296   //
1297   // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
1298   //
1299   //    NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
1300   //    The number of rows and COLUMNS can be different.
1301   //
1302   //    ROWS       ==  2**M + 1
1303   //    COLUMNS    ==  2**N + 1
1304   //    PHISLICES  ==  Arbitrary but greater than 3
1305   //
1306   //    DeltaPhi in Radians
1307   //
1308   //    SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
1309   //             = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
1310   //
1311   // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
1312
1313   const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
1314
1315   const Float_t  gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
1316   const Float_t  gridSizePhi =  deltaphi ;
1317   const Float_t  gridSizeZ   =  fgkTPCZ0 / (columns-1) ;
1318   const Float_t  ratioPhi    =  gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
1319   const Float_t  ratioZ      =  gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
1320
1321   TMatrixD arrayE(rows,columns) ;
1322
1323   // Check that the number of rows and columns is suitable for a binary expansion
1324   if ( !IsPowerOfTwo((rows-1))    ) {
1325     AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
1326     return; }
1327   if ( !IsPowerOfTwo((columns-1)) ) {
1328     AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
1329     return; }
1330   if ( phislices <= 3   )  {
1331     AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
1332     return; }
1333   if  ( phislices > 1000 ) {
1334     AliError("Poisson3D  phislices > 1000 is not allowed (nor wise) ");
1335     return; }
1336
1337   // Solve Poisson's equation in cylindrical coordinates by relaxation technique
1338   // Allow for different size grid spacing in R and Z directions
1339   // Use a binary expansion of the matrix to speed up the solution of the problem
1340
1341   Int_t loops, mplus, mminus, signplus, signminus  ;
1342   Int_t ione = (rows-1)/4 ;
1343   Int_t jone = (columns-1)/4 ;
1344   loops = TMath::Max(ione, jone) ;      // Calculate the number of loops for the binary expansion
1345   loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ;  // Solve for N in 2**N
1346
1347   TMatrixD* arrayofSumChargeDensities[1000] ;    // Create temporary arrays to store low resolution charge arrays
1348
1349   for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
1350
1351   for ( Int_t count = 0 ; count < loops ; count++ ) {      // START the master loop and do the binary expansion
1352
1353     Float_t  tempgridSizeR   =  gridSizeR  * ione ;
1354     Float_t  tempratioPhi    =  ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
1355     Float_t  tempratioZ      =  ratioZ   * ione * ione / ( jone * jone ) ;
1356
1357     std::vector<float> coef1(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1358     std::vector<float> coef2(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1359     std::vector<float> coef3(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1360     std::vector<float> coef4(rows) ;  // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1361
1362     for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
1363       Float_t radius = fgkIFCRadius + i*gridSizeR ;
1364       coef1[i] = 1.0 + tempgridSizeR/(2*radius);
1365       coef2[i] = 1.0 - tempgridSizeR/(2*radius);
1366       coef3[i] = tempratioPhi/(radius*radius);
1367       coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
1368     }
1369
1370     for ( Int_t m = 0 ; m < phislices ; m++ ) {
1371       TMatrixD &chargeDensity    = *arrayofChargeDensities[m] ;
1372       TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
1373       for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
1374         Float_t radius = fgkIFCRadius + i*gridSizeR ;
1375         for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
1376           if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1377           else {           // Add up all enclosed charge density contributions within 1/2 unit in all directions
1378             Float_t weight = 0.0 ;
1379             Float_t sum    = 0.0 ;
1380             sumChargeDensity(i,j) = 0.0 ;
1381             for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
1382               for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
1383                 if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
1384                 else
1385                   weight = 1.0 ;
1386                 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1387                 sum += weight*radius ;
1388               }
1389             }
1390             sumChargeDensity(i,j) /= sum ;
1391           }
1392           sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
1393         }
1394       }
1395     }
1396
1397     for ( Int_t k = 1 ; k <= iterations; k++ ) {
1398
1399       // over-relaxation index, >= 1 but < 2
1400       Float_t overRelax   = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1401       Float_t overRelaxM1 = overRelax - 1.0 ;
1402
1403       std::vector<float> overRelaxcoef4(rows) ;  // Do this the standard C++ way to avoid gcc extensions
1404       std::vector<float> overRelaxcoef5(rows) ;  // Do this the standard C++ way to avoid gcc extensions
1405
1406       for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1407         overRelaxcoef4[i] = overRelax * coef4[i] ;
1408         overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
1409       }
1410
1411       for ( Int_t m = 0 ; m < phislices ; m++ ) {
1412
1413         mplus  = m + 1;   signplus  = 1 ;
1414         mminus = m - 1 ;  signminus = 1 ;
1415         if (symmetry==1) {  // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1416           if ( mplus  > phislices-1 ) mplus  = phislices - 2 ;
1417           if ( mminus < 0 )           mminus = 1 ;
1418         }
1419         else if (symmetry==-1) {   // Anti-symmetry in phi
1420           if ( mplus  > phislices-1 ) { mplus  = phislices - 2 ; signplus  = -1 ; }
1421           if ( mminus < 0 )           { mminus = 1 ;             signminus = -1 ; }
1422         }
1423                 else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1424           if ( mplus  > phislices-1 ) mplus  = m + 1 - phislices ;
1425           if ( mminus < 0 )           mminus = m - 1 + phislices ;
1426         }
1427         TMatrixD& arrayV    =  *arrayofArrayV[m] ;
1428         TMatrixD& arrayVP   =  *arrayofArrayV[mplus] ;
1429         TMatrixD& arrayVM   =  *arrayofArrayV[mminus] ;
1430         TMatrixD& sumChargeDensity =  *arrayofSumChargeDensities[m] ;
1431
1432         for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
1433           for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1434
1435             arrayV(i,j) = (   coef2[i]          *   arrayV(i-ione,j)
1436                             + tempratioZ        * ( arrayV(i,j-jone)  +  arrayV(i,j+jone) )
1437                             - overRelaxcoef5[i] *   arrayV(i,j)
1438                             + coef1[i]          *   arrayV(i+ione,j)
1439                             + coef3[i]          * ( signplus*arrayVP(i,j)       +  signminus*arrayVM(i,j) )
1440                             + sumChargeDensity(i,j)
1441                           ) * overRelaxcoef4[i] ;
1442             // Note: over-relax the solution at each step.  This speeds up the convergance.
1443
1444           }
1445         }
1446
1447         if ( k == iterations ) {   // After full solution is achieved, copy low resolution solution into higher res array
1448           for ( Int_t i = ione ; i < rows-1 ; i+=ione )  {
1449             for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1450
1451               if ( ione > 1 ) {
1452                 arrayV(i+ione/2,j)                    =  ( arrayV(i+ione,j) + arrayV(i,j)     ) / 2 ;
1453                 if ( i == ione )  arrayV(i-ione/2,j) =  ( arrayV(0,j)       + arrayV(ione,j) ) / 2 ;
1454               }
1455               if ( jone > 1 ) {
1456                 arrayV(i,j+jone/2)                    =  ( arrayV(i,j+jone) + arrayV(i,j) )     / 2 ;
1457                 if ( j == jone )  arrayV(i,j-jone/2) =  ( arrayV(i,0)       + arrayV(i,jone) ) / 2 ;
1458               }
1459               if ( ione > 1 && jone > 1 ) {
1460                 arrayV(i+ione/2,j+jone/2) =  ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
1461                 if ( i == ione ) arrayV(i-ione/2,j-jone/2) =   ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
1462                 if ( j == jone ) arrayV(i-ione/2,j-jone/2) =   ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
1463                 // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
1464               }
1465             }
1466           }
1467         }
1468
1469       }
1470     }
1471
1472     ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
1473     jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
1474
1475   }
1476
1477   //Differentiate V(r) and solve for E(r) using special equations for the first and last row
1478   //Integrate E(r)/E(z) from point of origin to pad plane
1479
1480   for ( Int_t m = 0 ; m < phislices ; m++ ) {
1481     TMatrixD& arrayV    =  *arrayofArrayV[m] ;
1482     TMatrixD& eroverEz  =  *arrayofEroverEz[m] ;
1483
1484     for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {  // Count backwards to facilitate integration over Z
1485
1486       // Differentiate in R
1487       for ( Int_t i = 1 ; i < rows-1 ; i++ )  arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1488       arrayE(0,j)      =  -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1489       arrayE(rows-1,j) =  -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1490       // Integrate over Z
1491       for ( Int_t i = 0 ; i < rows ; i++ ) {
1492         Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
1493         eroverEz(i,j) = 0.0 ;
1494         for ( Int_t k = j ; k < columns ; k++ ) {
1495
1496           eroverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1497           if ( index != 4 )  index = 4; else index = 2 ;
1498         }
1499         if ( index == 4 ) eroverEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1500         if ( index == 2 ) eroverEz(i,j)  +=
1501           (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1502         if ( j == columns-2 ) eroverEz(i,j) =
1503           (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1504         if ( j == columns-1 ) eroverEz(i,j) =  0.0 ;
1505       }
1506     }
1507     // if ( m == 0 ) { TCanvas*  c1 =  new TCanvas("erOverEz","erOverEz",50,50,840,600) ;  c1 -> cd() ;
1508     // eroverEz.Draw("surf") ; } // JT test
1509   }
1510
1511   //Differentiate V(r) and solve for E(phi)
1512   //Integrate E(phi)/E(z) from point of origin to pad plane
1513
1514   for ( Int_t m = 0 ; m < phislices ; m++ ) {
1515
1516     mplus  = m + 1;   signplus  = 1 ;
1517     mminus = m - 1 ;  signminus = 1 ;
1518     if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1519       if ( mplus  > phislices-1 ) mplus  = phislices - 2 ;
1520       if ( mminus < 0 )           mminus = 1 ;
1521     }
1522     else if (symmetry==-1) {       // Anti-symmetry in phi
1523       if ( mplus  > phislices-1 ) { mplus  = phislices - 2 ;  signplus  = -1 ; }
1524       if ( mminus < 0 )           { mminus = 1 ;                    signminus = -1 ; }
1525     }
1526     else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
1527       if ( mplus  > phislices-1 ) mplus  = m + 1 - phislices ;
1528       if ( mminus < 0 )           mminus = m - 1 + phislices ;
1529     }
1530     TMatrixD &arrayVP     =  *arrayofArrayV[mplus] ;
1531     TMatrixD &arrayVM     =  *arrayofArrayV[mminus] ;
1532     TMatrixD &ePhioverEz  =  *arrayofEPhioverEz[m] ;
1533     for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1534       // Differentiate in Phi
1535       for ( Int_t i = 0 ; i < rows ; i++ ) {
1536         Float_t radius = fgkIFCRadius + i*gridSizeR ;
1537         arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
1538       }
1539       // Integrate over Z
1540       for ( Int_t i = 0 ; i < rows ; i++ ) {
1541         Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
1542         ePhioverEz(i,j) = 0.0 ;
1543         for ( Int_t k = j ; k < columns ; k++ ) {
1544
1545           ePhioverEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1546           if ( index != 4 )  index = 4; else index = 2 ;
1547         }
1548         if ( index == 4 ) ePhioverEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1549         if ( index == 2 ) ePhioverEz(i,j)  +=
1550           (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1551         if ( j == columns-2 ) ePhioverEz(i,j) =
1552           (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1553         if ( j == columns-1 ) ePhioverEz(i,j) =  0.0 ;
1554       }
1555     }
1556     // if ( m == 5 ) { TCanvas* c2 =  new TCanvas("arrayE","arrayE",50,50,840,600) ;  c2 -> cd() ;
1557     // arrayE.Draw("surf") ; } // JT test
1558   }
1559
1560
1561   // Differentiate V(r) and solve for E(z) using special equations for the first and last row
1562   // Integrate (E(z)-Ezstd) from point of origin to pad plane
1563
1564   for ( Int_t m = 0 ; m < phislices ; m++ ) {
1565     TMatrixD& arrayV   =  *arrayofArrayV[m] ;
1566     TMatrixD& deltaEz  =  *arrayofDeltaEz[m] ;
1567
1568     // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1569     for ( Int_t i = 0 ; i < rows ; i++) {
1570       for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1571       arrayE(i,0)         =  -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1572       arrayE(i,columns-1) =  -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1573     }
1574
1575     for ( Int_t j = columns-1 ; j >= 0 ; j-- ) {  // Count backwards to facilitate integration over Z
1576       // Integrate over Z
1577       for ( Int_t i = 0 ; i < rows ; i++ ) {
1578         Int_t index = 1 ;   // Simpsons rule if N=odd.  If N!=odd then add extra point by trapezoidal rule.
1579         deltaEz(i,j) = 0.0 ;
1580         for ( Int_t k = j ; k < columns ; k++ ) {
1581           deltaEz(i,j)  +=  index*(gridSizeZ/3.0)*arrayE(i,k) ;
1582           if ( index != 4 )  index = 4; else index = 2 ;
1583         }
1584         if ( index == 4 ) deltaEz(i,j)  -=  (gridSizeZ/3.0)*arrayE(i,columns-1) ;
1585         if ( index == 2 ) deltaEz(i,j)  +=
1586           (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
1587         if ( j == columns-2 ) deltaEz(i,j) =
1588           (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
1589         if ( j == columns-1 ) deltaEz(i,j) =  0.0 ;
1590       }
1591     }
1592     // if ( m == 0 ) { TCanvas*  c1 =  new TCanvas("erOverEz","erOverEz",50,50,840,600) ;  c1 -> cd() ;
1593     // eroverEz.Draw("surf") ; } // JT test
1594
1595     // calculate z distortion from the integrated Delta Ez residuals
1596     // and include the aquivalence (Volt to cm) of the ROC shift !!
1597
1598     for ( Int_t j = 0 ; j < columns ; j++ )  {
1599       for ( Int_t i = 0 ; i < rows ; i++ ) {
1600
1601         // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1602         deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
1603
1604         // ROC Potential in cm aquivalent
1605         Double_t dzROCShift =  arrayV(i, columns -1)/ezField;
1606         if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift;  // add the ROC misaligment
1607
1608       }
1609     }
1610
1611   } // end loop over phi
1612
1613
1614
1615   for ( Int_t k = 0 ; k < phislices ; k++ )
1616     {
1617       arrayofSumChargeDensities[k]->Delete() ;
1618     }
1619
1620
1621
1622   arrayE.Clear();
1623 }
1624
1625
1626 Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1627   //
1628   // Helperfunction: Check if integer is a power of 2
1629   //
1630   Int_t j = 0;
1631   while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
1632   if ( j == 1 ) return(1) ;  // True
1633   return(0) ;                // False
1634 }
1635
1636
1637 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
1638   //
1639   // Fit the track parameters - without and with distortion
1640   // 1. Space points in the TPC are simulated along the trajectory
1641   // 2. Space points distorted
1642   // 3. Fits  the non distorted and distroted track to the reference plane at refX
1643   // 4. For visualization and debugging  purposes the space points and tracks can be stored  in the tree - using the TTreeSRedirector functionality
1644   //
1645   // trackIn   - input track parameters
1646   // refX     - reference X to fit the track
1647   // dir      - direction - out=1 or in=-1
1648   // pcstream -  debug streamer to check the results
1649   //
1650   // see AliExternalTrackParam.h documentation:
1651   // track1.fP[0] - local y (rphi)
1652   // track1.fP[1] - z
1653   // track1.fP[2] - sinus of local inclination angle
1654   // track1.fP[3] - tangent of deep angle
1655   // track1.fP[4] - 1/pt
1656
1657   AliTPCROC * roc = AliTPCROC::Instance();
1658   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
1659   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
1660   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1661   const Double_t kMaxSnp = 0.85;
1662   const Double_t kSigmaY=0.1;
1663   const Double_t kSigmaZ=0.1;
1664   const Double_t kMaxR=500;
1665   const Double_t kMaxZ=500;
1666
1667   const Double_t kMaxZ0=220;
1668   const Double_t kZcut=3;
1669   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1670   Int_t npoints1=0;
1671   Int_t npoints2=0;
1672
1673   AliExternalTrackParam  track(trackIn); //
1674   // generate points
1675   AliTrackPointArray pointArray0(npoints0);
1676   AliTrackPointArray pointArray1(npoints0);
1677   Double_t xyz[3];
1678   if (!AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,5,kTRUE,kMaxSnp)) return 0;
1679   //
1680   // simulate the track
1681   Int_t npoints=0;
1682   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
1683   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
1684     if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) return 0;
1685     track.GetXYZ(xyz);
1686     xyz[0]+=gRandom->Gaus(0,0.000005);
1687     xyz[1]+=gRandom->Gaus(0,0.000005);
1688     xyz[2]+=gRandom->Gaus(0,0.000005);
1689     if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
1690     if (TMath::Abs(track.GetX())<kRTPC0) continue;
1691     if (TMath::Abs(track.GetX())>kRTPC1) continue;
1692     AliTrackPoint pIn0;                               // space point
1693     AliTrackPoint pIn1;
1694     Int_t sector= (xyz[2]>0)? 0:18;
1695     pointArray0.GetPoint(pIn0,npoints);
1696     pointArray1.GetPoint(pIn1,npoints);
1697     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
1698     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
1699     DistortPoint(distPoint, sector);
1700     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
1701     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
1702     //
1703     track.Rotate(alpha);
1704     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
1705     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
1706     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
1707     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
1708     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
1709     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
1710     pointArray0.AddPoint(npoints, &pIn0);
1711     pointArray1.AddPoint(npoints, &pIn1);
1712     npoints++;
1713     if (npoints>=npoints0) break;
1714   }
1715   if (npoints<npoints0/4.) return 0;
1716   //
1717   // refit track
1718   //
1719   AliExternalTrackParam *track0=0;
1720   AliExternalTrackParam *track1=0;
1721   AliTrackPoint   point1,point2,point3;
1722   if (dir==1) {  //make seed inner
1723     pointArray0.GetPoint(point1,1);
1724     pointArray0.GetPoint(point2,11);
1725     pointArray0.GetPoint(point3,21);
1726   }
1727   if (dir==-1){ //make seed outer
1728     pointArray0.GetPoint(point1,npoints-21);
1729     pointArray0.GetPoint(point2,npoints-11);
1730     pointArray0.GetPoint(point3,npoints-1);
1731   }
1732   if ((TMath::Abs(point1.GetX()-point3.GetX())+TMath::Abs(point1.GetY()-point3.GetY()))<10){
1733     printf("fit points not properly initialized\n");
1734     return 0;
1735   }
1736   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
1737   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
1738   track0->ResetCovariance(10);
1739   track1->ResetCovariance(10);
1740   if (TMath::Abs(AliTrackerBase::GetBz())<0.01){
1741     ((Double_t*)track0->GetParameter())[4]=  trackIn.GetParameter()[4];
1742     ((Double_t*)track1->GetParameter())[4]=  trackIn.GetParameter()[4];
1743   }
1744   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
1745     Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
1746     //
1747     AliTrackPoint pIn0;
1748     AliTrackPoint pIn1;
1749     pointArray0.GetPoint(pIn0,ipoint);
1750     pointArray1.GetPoint(pIn1,ipoint);
1751     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
1752     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
1753     if (TMath::Abs(prot0.GetX())<kRTPC0) continue;
1754     if (TMath::Abs(prot0.GetX())>kRTPC1) continue;
1755     //
1756     if (!AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
1757     if (!AliTrackerBase::PropagateTrackTo(track1,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
1758     if (TMath::Abs(track0->GetZ())>kMaxZ) break;
1759     if (TMath::Abs(track0->GetX())>kMaxR) break;
1760     if (TMath::Abs(track1->GetZ())>kMaxZ) break;
1761     if (TMath::Abs(track1->GetX())>kMaxR) break;
1762     if (dir>0 && track1->GetX()>refX) continue;
1763     if (dir<0 && track1->GetX()<refX) continue;
1764     if (TMath::Abs(track1->GetZ())<kZcut)continue;
1765     track.GetXYZ(xyz);  // distorted track also propagated to the same reference radius
1766     //
1767     Double_t pointPos[2]={0,0};
1768     Double_t pointCov[3]={0,0,0};
1769     pointPos[0]=prot0.GetY();//local y
1770     pointPos[1]=prot0.GetZ();//local z
1771     pointCov[0]=prot0.GetCov()[3];//simay^2
1772     pointCov[1]=prot0.GetCov()[4];//sigmayz
1773     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
1774     if (!track0->Update(pointPos,pointCov)) break;
1775     //
1776     Double_t deltaX=prot1.GetX()-prot0.GetX();   // delta X
1777     Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp()));  // deltaY due  delta X
1778     Double_t deltaZX=deltaX*track1->GetTgl();                           // deltaZ due  delta X
1779
1780     pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
1781     pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
1782     pointCov[0]=prot1.GetCov()[3];//simay^2
1783     pointCov[1]=prot1.GetCov()[4];//sigmayz
1784     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
1785     if (!track1->Update(pointPos,pointCov)) break;
1786     npoints1++;
1787     npoints2++;
1788   }
1789   if (npoints2<npoints/4.)  return 0;
1790   AliTrackerBase::PropagateTrackTo(track0,refX,kMass,5.,kTRUE,kMaxSnp);
1791   AliTrackerBase::PropagateTrackTo(track0,refX,kMass,1.,kTRUE,kMaxSnp);
1792   track1->Rotate(track0->GetAlpha());
1793   AliTrackerBase::PropagateTrackTo(track1,track0->GetX(),kMass,5.,kFALSE,kMaxSnp);
1794
1795   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
1796     "point0.="<<&pointArray0<<   //  points
1797     "point1.="<<&pointArray1<<   //  distorted points
1798     "trackIn.="<<&track<<       //  original track
1799     "track0.="<<track0<<         //  fitted track
1800     "track1.="<<track1<<         //  fitted distorted track
1801     "\n";
1802   new(&trackIn) AliExternalTrackParam(*track0);
1803   delete track0;
1804   return track1;
1805 }
1806
1807
1808
1809
1810
1811 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
1812   //
1813   // create the distortion tree on a mesh with granularity given by step
1814   // return the tree with distortions at given position
1815   // Map is created on the mesh with given step size
1816   //
1817   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
1818   Float_t xyz[3];
1819   for (Double_t x= -250; x<250; x+=step){
1820     for (Double_t y= -250; y<250; y+=step){
1821       Double_t r    = TMath::Sqrt(x*x+y*y);
1822       if (r<80) continue;
1823       if (r>250) continue;
1824       for (Double_t z= -250; z<250; z+=step){
1825         Int_t roc=(z>0)?0:18;
1826         xyz[0]=x;
1827         xyz[1]=y;
1828         xyz[2]=z;
1829         Double_t phi  = TMath::ATan2(y,x);
1830         DistortPoint(xyz,roc);
1831         Double_t r1    = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1832         Double_t phi1  = TMath::ATan2(xyz[1],xyz[0]);
1833         if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
1834         if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
1835         Double_t dx = xyz[0]-x;
1836         Double_t dy = xyz[1]-y;
1837         Double_t dz = xyz[2]-z;
1838         Double_t dr=r1-r;
1839         Double_t drphi=(phi1-phi)*r;
1840         (*pcstream)<<"distortion"<<
1841           "x="<<x<<           // original position
1842           "y="<<y<<
1843           "z="<<z<<
1844           "r="<<r<<
1845           "phi="<<phi<<
1846           "x1="<<xyz[0]<<      // distorted position
1847           "y1="<<xyz[1]<<
1848           "z1="<<xyz[2]<<
1849           "r1="<<r1<<
1850           "phi1="<<phi1<<
1851           //
1852           "dx="<<dx<<          // delta position
1853           "dy="<<dy<<
1854           "dz="<<dz<<
1855           "dr="<<dr<<
1856           "drphi="<<drphi<<
1857           "\n";
1858       }
1859     }
1860   }
1861   delete pcstream;
1862   TFile f(Form("correction%s.root",GetName()));
1863   TTree * tree = (TTree*)f.Get("distortion");
1864   TTree * tree2= tree->CopyTree("1");
1865   tree2->SetName(Form("dist%s",GetName()));
1866   tree2->SetDirectory(0);
1867   delete tree;
1868   return tree2;
1869 }
1870
1871
1872
1873
1874 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1875   //
1876   // Make a fit tree:
1877   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1878   // calculates partial distortions
1879   // Partial distortion is stored in the resulting tree
1880   // Output is storred in the file distortion_<dettype>_<partype>.root
1881   // Partial  distortion is stored with the name given by correction name
1882   //
1883   //
1884   // Parameters of function:
1885   // input     - input tree
1886   // dtype     - distortion type 0 - ITSTPC,  1 -TPCTRD, 2 - TPCvertex , 3 - TPC-TOF,  4 - TPCTPC track crossing
1887   // ppype     - parameter type
1888   // corrArray - array with partial corrections
1889   // step      - skipe entries  - if 1 all entries processed - it is slow
1890   // debug     0 if debug on also space points dumped - it is slow
1891
1892   const Double_t kMaxSnp = 0.85;
1893   const Double_t kcutSnp=0.25;
1894   const Double_t kcutTheta=1.;
1895   const Double_t kRadiusTPC=85;
1896   //  AliTPCROC *tpcRoc =AliTPCROC::Instance();
1897   //
1898   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1899   //  const Double_t kB2C=-0.299792458e-3;
1900   const Int_t kMinEntries=20;
1901   Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
1902   Float_t refX;
1903   Int_t run;
1904   tinput->SetBranchAddress("run",&run);
1905   tinput->SetBranchAddress("theta",&theta);
1906   tinput->SetBranchAddress("phi", &phi);
1907   tinput->SetBranchAddress("snp",&snp);
1908   tinput->SetBranchAddress("mean",&mean);
1909   tinput->SetBranchAddress("rms",&rms);
1910   tinput->SetBranchAddress("entries",&entries);
1911   tinput->SetBranchAddress("sector",&sector);
1912   tinput->SetBranchAddress("dsec",&dsec);
1913   tinput->SetBranchAddress("refX",&refX);
1914   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
1915   //
1916   Int_t nentries=tinput->GetEntries();
1917   Int_t ncorr=corrArray->GetEntries();
1918   Double_t corrections[100]={0}; //
1919   Double_t tPar[5];
1920   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1921   Int_t dir=0;
1922   if (dtype==5 || dtype==6) dtype=4;
1923   if (dtype==0) { dir=-1;}
1924   if (dtype==1) { dir=1;}
1925   if (dtype==2) { dir=-1;}
1926   if (dtype==3) { dir=1;}
1927   if (dtype==4) { dir=-1;}
1928   //
1929   for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1930     tinput->GetEntry(ientry);
1931     if (TMath::Abs(snp)>kMaxSnp) continue;
1932     tPar[0]=0;
1933     tPar[1]=theta*refX;
1934     if (dtype==2)  tPar[1]=theta*kRadiusTPC;
1935     tPar[2]=snp;
1936     tPar[3]=theta;
1937     tPar[4]=(gRandom->Rndm()-0.5)*0.02;  // should be calculated - non equal to 0
1938     if (dtype==4){
1939       // tracks crossing CE
1940       tPar[1]=0;   // track at the CE
1941       //if (TMath::Abs(theta) <0.05) continue;  // deep cross
1942     }
1943
1944     if (TMath::Abs(snp) >kcutSnp) continue;
1945     if (TMath::Abs(theta) >kcutTheta) continue;
1946     printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1947     Double_t bz=AliTrackerBase::GetBz();
1948     if (dtype !=4) { //exclude TPC  - for TPC mainly non primary tracks
1949       if (dtype!=2 && TMath::Abs(bz)>0.1 )  tPar[4]=snp/(refX*bz*kB2C*2);
1950
1951       if (dtype==2 && TMath::Abs(bz)>0.1 )  {
1952         tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);//
1953         // snp at the TPC inner radius in case the vertex match used
1954       }
1955     }
1956     //
1957     tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
1958     AliExternalTrackParam track(refX,phi,tPar,cov);
1959     Double_t xyz[3];
1960     track.GetXYZ(xyz);
1961     Int_t id=0;
1962     Double_t pt=1./tPar[4];
1963     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1964     //if (ptype==4 &&bz<0) mean*=-1;  // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used
1965     Double_t refXD=refX;
1966     (*pcstream)<<"fit"<<
1967       "run="<<run<<       // run number
1968       "bz="<<bz<<         // magnetic filed used
1969       "dtype="<<dtype<<   // detector match type
1970       "ptype="<<ptype<<   // parameter type
1971       "theta="<<theta<<   // theta
1972       "phi="<<phi<<       // phi
1973       "snp="<<snp<<       // snp
1974       "mean="<<mean<<     // mean dist value
1975       "rms="<<rms<<       // rms
1976       "sector="<<sector<<
1977       "dsec="<<dsec<<
1978       "refX="<<refXD<<         // referece X as double
1979       "gx="<<xyz[0]<<         // global position at reference
1980       "gy="<<xyz[1]<<         // global position at reference
1981       "gz="<<xyz[2]<<         // global position at reference
1982       "dRrec="<<dRrec<<      // delta Radius in reconstruction
1983       "pt="<<pt<<            // pt
1984       "id="<<id<<             // track id
1985       "entries="<<entries;// number of entries in bin
1986     //
1987     Bool_t isOK=kTRUE;
1988     if (entries<kMinEntries) isOK=kFALSE;
1989     //
1990     if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1991       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1992       corrections[icorr]=0;
1993       if (entries>kMinEntries){
1994         AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1995         AliExternalTrackParam *trackOut = 0;
1996         if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1997         if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1998         if (dtype==0) {dir= -1;}
1999         if (dtype==1) {dir=  1;}
2000         if (dtype==2) {dir= -1;}
2001         if (dtype==3) {dir=  1;}
2002         //
2003         if (trackOut){
2004           if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
2005           if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE;
2006           if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
2007           //      trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
2008           //
2009           corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
2010           delete trackOut;
2011         }else{
2012           corrections[icorr]=0;
2013           isOK=kFALSE;
2014         }
2015         //if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature - commented out
2016       }
2017       (*pcstream)<<"fit"<<
2018         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
2019     }
2020
2021     if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
2022       //
2023       // special case of the TPC tracks crossing the CE
2024       //
2025       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
2026       corrections[icorr]=0;
2027       if (entries>kMinEntries){
2028         AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
2029         AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet
2030         AliExternalTrackParam *trackOut0 = 0;
2031         AliExternalTrackParam *trackOut1 = 0;
2032         //
2033         if (debug)  trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
2034         if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
2035         if (debug)  trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
2036         if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
2037         //
2038         if (trackOut0 && trackOut1){
2039           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp))  isOK=kFALSE;
2040           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2041           if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2042           if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
2043           //
2044           if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
2045           if (!trackIn1.Rotate(trackIn0.GetAlpha()))  isOK=kFALSE;
2046           if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2047           if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE;
2048           if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
2049           //
2050           corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
2051           corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
2052           if (isOK)
2053             if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
2054                 (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
2055                 (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
2056                 (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
2057                 (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
2058                 (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
2059                 ){
2060               isOK=kFALSE;
2061             }
2062           delete trackOut0;
2063           delete trackOut1;
2064         }else{
2065           corrections[icorr]=0;
2066           isOK=kFALSE;
2067         }
2068         //
2069         //if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature - commented out no in lookup
2070       }
2071       (*pcstream)<<"fit"<<
2072         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
2073     }
2074     //
2075     (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
2076   }
2077
2078
2079   delete pcstream;
2080 }
2081
2082
2083
2084 void AliTPCCorrection::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
2085   //
2086   // Make a fit tree:
2087   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
2088   // calculates partial distortions
2089   // Partial distortion is stored in the resulting tree
2090   // Output is storred in the file distortion_<dettype>_<partype>.root
2091   // Partial  distortion is stored with the name given by correction name
2092   //
2093   //
2094   // Parameters of function:
2095   // input     - input tree
2096   // dtype     - distortion type 10 - IROC-OROC
2097   // ppype     - parameter type
2098   // corrArray - array with partial corrections
2099   // step      - skipe entries  - if 1 all entries processed - it is slow
2100   // debug     0 if debug on also space points dumped - it is slow
2101
2102   const Double_t kMaxSnp = 0.8;
2103   const Int_t kMinEntries=200;
2104   //  AliTPCROC *tpcRoc =AliTPCROC::Instance();
2105   //
2106   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2107   //  const Double_t kB2C=-0.299792458e-3;
2108   Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
2109   Int_t isec1, isec0;
2110   Double_t refXD;
2111   Float_t refX;
2112   Int_t run;
2113   tinput->SetBranchAddress("run",&run);
2114   tinput->SetBranchAddress("theta",&theta);
2115   tinput->SetBranchAddress("phi", &phi);
2116   tinput->SetBranchAddress("snp",&snp);
2117   tinput->SetBranchAddress("mean",&mean);
2118   tinput->SetBranchAddress("rms",&rms);
2119   tinput->SetBranchAddress("entries",&entries);
2120   tinput->SetBranchAddress("sector",&sector);
2121   tinput->SetBranchAddress("dsec",&dsec);
2122   tinput->SetBranchAddress("refX",&refXD);
2123   tinput->SetBranchAddress("z",&globalZ);
2124   tinput->SetBranchAddress("isec0",&isec0);
2125   tinput->SetBranchAddress("isec1",&isec1);
2126   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
2127   //
2128   Int_t nentries=tinput->GetEntries();
2129   Int_t ncorr=corrArray->GetEntries();
2130   Double_t corrections[100]={0}; //
2131   Double_t tPar[5];
2132   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
2133   Int_t dir=0;
2134   //
2135   for (Int_t ientry=offset; ientry<nentries; ientry+=step){
2136     tinput->GetEntry(ientry);
2137     refX=refXD;
2138     Int_t id=-1;
2139     if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1;  // IROC-OROC - opposite side
2140     if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2;  // IROC-OROC - same side
2141     if (dtype==10  && id==-1) continue;
2142     //
2143     dir=-1;
2144     tPar[0]=0;
2145     tPar[1]=globalZ;
2146     tPar[2]=snp;
2147     tPar[3]=theta;
2148     tPar[4]=(gRandom->Rndm()-0.1)*0.2;  //
2149     Double_t pt=1./tPar[4];
2150     //
2151     printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
2152     Double_t bz=AliTrackerBase::GetBz();
2153     AliExternalTrackParam track(refX,phi,tPar,cov);
2154     Double_t xyz[3],xyzIn[3],xyzOut[3];
2155     track.GetXYZ(xyz);
2156     track.GetXYZAt(85,bz,xyzIn);
2157     track.GetXYZAt(245,bz,xyzOut);
2158     Double_t phiIn  = TMath::ATan2(xyzIn[1],xyzIn[0]);
2159     Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
2160     Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
2161     Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
2162     Int_t sectorIn  = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
2163     Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
2164     //
2165     Bool_t isOK=kTRUE;
2166     if (sectorIn!=sectorOut) isOK=kFALSE;  // requironment - cluster in the same sector
2167     if (sectorIn!=sectorRef) isOK=kFALSE;  // requironment - cluster in the same sector
2168     if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE;  // requironment - minimal amount of tracks in bin
2169     // Do downscale
2170     if (TMath::Abs(theta)>1) isOK=kFALSE;
2171     //
2172     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
2173     //
2174     (*pcstream)<<"fit"<<
2175       "run="<<run<<       //run
2176       "bz="<<bz<<         // magnetic filed used
2177       "dtype="<<dtype<<   // detector match type
2178       "ptype="<<ptype<<   // parameter type
2179       "theta="<<theta<<   // theta
2180       "phi="<<phi<<       // phi
2181       "snp="<<snp<<       // snp
2182       "mean="<<mean<<     // mean dist value
2183       "rms="<<rms<<       // rms
2184       "sector="<<sector<<
2185       "dsec="<<dsec<<
2186       "refX="<<refXD<<         // referece X
2187       "gx="<<xyz[0]<<         // global position at reference
2188       "gy="<<xyz[1]<<         // global position at reference
2189       "gz="<<xyz[2]<<         // global position at reference
2190       "dRrec="<<dRrec<<      // delta Radius in reconstruction
2191       "pt="<<pt<<      //pt
2192       "id="<<id<<             // track id
2193       "entries="<<entries;// number of entries in bin
2194     //
2195     AliExternalTrackParam *trackOut0 = 0;
2196     AliExternalTrackParam *trackOut1 = 0;
2197     AliExternalTrackParam *ptrackIn0 = 0;
2198     AliExternalTrackParam *ptrackIn1 = 0;
2199
2200     for (Int_t icorr=0; icorr<ncorr; icorr++) {
2201       //
2202       // special case of the TPC tracks crossing the CE
2203       //
2204       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
2205       corrections[icorr]=0;
2206       if (entries>kMinEntries &&isOK){
2207         AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
2208         AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
2209         ptrackIn1=&trackIn0;
2210         ptrackIn0=&trackIn1;
2211         //
2212         if (debug)  trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
2213         if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
2214         if (debug)  trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
2215         if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
2216         //
2217         if (trackOut0 && trackOut1){
2218           //
2219           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp))  isOK=kFALSE;
2220           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2221           // rotate all tracks to the same frame
2222           if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2223           if (!trackIn1.Rotate(trackIn0.GetAlpha()))  isOK=kFALSE;
2224           if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2225           //
2226           if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2227           if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2228           if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2229           //
2230           corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
2231           corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
2232           (*pcstream)<<"fitDebug"<< // just to debug the correction
2233             "mean="<<mean<<
2234             "pIn0.="<<ptrackIn0<<
2235             "pIn1.="<<ptrackIn1<<
2236             "pOut0.="<<trackOut0<<
2237             "pOut1.="<<trackOut1<<
2238             "refX="<<refXD<<
2239             "\n";
2240           delete trackOut0;
2241           delete trackOut1;
2242         }else{
2243           corrections[icorr]=0;
2244           isOK=kFALSE;
2245         }
2246       }
2247       (*pcstream)<<"fit"<<
2248         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
2249     }
2250     //
2251     (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
2252   }
2253   delete pcstream;
2254 }
2255
2256
2257
2258 void AliTPCCorrection::MakeLaserDistortionTreeOld(TTree* tree, TObjArray *corrArray, Int_t itype){
2259   //
2260   // Make a laser fit tree for global minimization
2261   //
2262   const Double_t cutErrY=0.1;
2263   const Double_t cutErrZ=0.1;
2264   const Double_t kEpsilon=0.00000001;
2265   const Double_t kMaxDist=1.;  // max distance - space correction
2266   const Double_t kMaxRMS=0.05;  // max distance -between point and local mean
2267   TVectorD *vecdY=0;
2268   TVectorD *vecdZ=0;
2269   TVectorD *veceY=0;
2270   TVectorD *veceZ=0;
2271   AliTPCLaserTrack *ltr=0;
2272   AliTPCLaserTrack::LoadTracks();
2273   tree->SetBranchAddress("dY.",&vecdY);
2274   tree->SetBranchAddress("dZ.",&vecdZ);
2275   tree->SetBranchAddress("eY.",&veceY);
2276   tree->SetBranchAddress("eZ.",&veceZ);
2277   tree->SetBranchAddress("LTr.",&ltr);
2278   Int_t entries= tree->GetEntries();
2279   TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
2280   Double_t bz=AliTrackerBase::GetBz();
2281   //
2282
2283   for (Int_t ientry=0; ientry<entries; ientry++){
2284     tree->GetEntry(ientry);
2285     if (!ltr->GetVecGX()){
2286       ltr->UpdatePoints();
2287     }
2288     TVectorD * delta= (itype==0)? vecdY:vecdZ;
2289     TVectorD * err= (itype==0)? veceY:veceZ;
2290     TLinearFitter  fitter(2,"pol1");
2291     for (Int_t iter=0; iter<2; iter++){
2292       Double_t kfit0=0, kfit1=0;
2293       Int_t npoints=fitter.GetNpoints();
2294       if (npoints>80){
2295         fitter.Eval();
2296         kfit0=fitter.GetParameter(0);
2297         kfit1=fitter.GetParameter(1);
2298       }
2299       for (Int_t irow=0; irow<159; irow++){
2300         Bool_t isOK=kTRUE;
2301         Int_t isOKF=0;
2302         Int_t nentries = 1000;
2303         if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
2304         if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
2305         Int_t dtype=5;
2306         Double_t array[10];
2307         Int_t first3=TMath::Max(irow-3,0);
2308         Int_t last3 =TMath::Min(irow+3,159);
2309         Int_t counter=0;
2310         if ((*ltr->GetVecSec())[irow]>=0 && err) {
2311           for (Int_t jrow=first3; jrow<=last3; jrow++){
2312             if ((*ltr->GetVecSec())[irow]!= (*ltr->GetVecSec())[jrow]) continue;
2313             if ((*err)[jrow]<kEpsilon) continue;
2314             array[counter]=(*delta)[jrow];
2315             counter++;
2316           }
2317         }
2318         Double_t rms3  = 0;
2319         Double_t mean3 = 0;
2320         if (counter>2){
2321           rms3  = TMath::RMS(counter,array);
2322           mean3  = TMath::Mean(counter,array);
2323         }else{
2324           isOK=kFALSE;
2325         }
2326         Double_t phi   =(*ltr->GetVecPhi())[irow];
2327         Double_t theta =ltr->GetTgl();
2328         Double_t mean=delta->GetMatrixArray()[irow];
2329         Double_t gx=0,gy=0,gz=0;
2330         Double_t snp = (*ltr->GetVecP2())[irow];
2331         Double_t dRrec=0;
2332         //      Double_t rms = err->GetMatrixArray()[irow];
2333         //
2334         gx = (*ltr->GetVecGX())[irow];
2335         gy = (*ltr->GetVecGY())[irow];
2336         gz = (*ltr->GetVecGZ())[irow];
2337         //
2338         // get delta R used in reconstruction
2339         AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
2340         AliTPCCorrection * correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
2341         //      const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
2342         //Double_t xyz0[3]={gx,gy,gz};
2343         Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
2344         Double_t fphi = TMath::ATan2(gy,gx);
2345         Double_t fsector = 9.*fphi/TMath::Pi();
2346         if (fsector<0) fsector+=18;
2347         Double_t dsec = fsector-Int_t(fsector)-0.5;
2348         Double_t refX=0;
2349         Int_t id= ltr->GetId();
2350         Double_t pt=0;
2351         //
2352         if (1 && oldR>1) {
2353           Float_t xyz1[3]={gx,gy,gz};
2354           Int_t sector=(gz>0)?0:18;
2355           correction->CorrectPoint(xyz1, sector);
2356           refX=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
2357           dRrec=oldR-refX;
2358         }
2359         if (TMath::Abs(rms3)>kMaxRMS) isOK=kFALSE;
2360         if (TMath::Abs(mean-mean3)>kMaxRMS) isOK=kFALSE;
2361         if (counter<4) isOK=kFALSE;
2362         if (npoints<90) isOK=kFALSE;
2363         if (isOK){
2364           fitter.AddPoint(&refX,mean);
2365         }
2366         Double_t deltaF=kfit0+kfit1*refX;
2367         if (iter==1){
2368           (*pcstream)<<"fitFull"<<  // dumpe also intermediate results
2369             "bz="<<bz<<         // magnetic filed used
2370             "dtype="<<dtype<<   // detector match type
2371             "ptype="<<itype<<   // parameter type
2372             "theta="<<theta<<   // theta
2373             "phi="<<phi<<       // phi
2374             "snp="<<snp<<       // snp
2375             "mean="<<mean3<<     // mean dist value
2376             "rms="<<rms3<<       // rms
2377             "deltaF="<<deltaF<<
2378             "npoints="<<npoints<<  //number of points
2379             "mean3="<<mean3<<     // mean dist value
2380             "rms3="<<rms3<<       // rms
2381             "counter="<<counter<<
2382             "sector="<<fsector<<
2383             "dsec="<<dsec<<
2384             //
2385             "refX="<<refX<<      // reference radius
2386             "gx="<<gx<<         // global position
2387             "gy="<<gy<<         // global position
2388             "gz="<<gz<<         // global position
2389             "dRrec="<<dRrec<<      // delta Radius in reconstruction
2390             "id="<<id<<     //bundle
2391             "entries="<<nentries<<// number of entries in bin
2392             "\n";
2393         }
2394         if (iter==1) (*pcstream)<<"fit"<<  // dump valus for fit
2395           "bz="<<bz<<         // magnetic filed used
2396           "dtype="<<dtype<<   // detector match type
2397           "ptype="<<itype<<   // parameter type
2398           "theta="<<theta<<   // theta
2399           "phi="<<phi<<       // phi
2400           "snp="<<snp<<       // snp
2401           "mean="<<mean3<<     // mean dist value
2402           "rms="<<rms3<<       // rms
2403           "sector="<<fsector<<
2404           "dsec="<<dsec<<
2405           //
2406           "refX="<<refX<<      // reference radius
2407           "gx="<<gx<<         // global position
2408           "gy="<<gy<<         // global position
2409           "gz="<<gz<<         // global position
2410           "dRrec="<<dRrec<<      // delta Radius in reconstruction
2411           "pt="<<pt<<           //pt
2412           "id="<<id<<     //bundle
2413           "entries="<<nentries;// number of entries in bin
2414         //
2415         //
2416         Double_t ky = TMath::Tan(TMath::ASin(snp));
2417         Int_t ncorr = corrArray->GetEntries();
2418         Double_t r0   = TMath::Sqrt(gx*gx+gy*gy);
2419         Double_t phi0 = TMath::ATan2(gy,gx);
2420         Double_t distortions[1000]={0};
2421         Double_t distortionsR[1000]={0};
2422         if (iter==1){
2423           for (Int_t icorr=0; icorr<ncorr; icorr++) {
2424             AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
2425             Float_t distPoint[3]={gx,gy,gz};
2426             Int_t sector= (gz>0)? 0:18;
2427             if (r0>80){
2428               corr->DistortPoint(distPoint, sector);
2429             }
2430             // Double_t value=distPoint[2]-gz;
2431             if (itype==0 && r0>1){
2432               Double_t r1   = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
2433               Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
2434               Double_t drphi= r0*(phi1-phi0);
2435               Double_t dr   = r1-r0;
2436               distortions[icorr]  = drphi-ky*dr;
2437               distortionsR[icorr] = dr;
2438             }
2439             if (TMath::Abs(distortions[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE; }
2440             if (TMath::Abs(distortionsR[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE;}
2441             (*pcstream)<<"fit"<<
2442               Form("%s=",corr->GetName())<<distortions[icorr];    // dump correction value
2443           }
2444           (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
2445         }
2446       }
2447     }
2448   }
2449   delete pcstream;
2450 }
2451
2452
2453
2454 void   AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type, Int_t integ){
2455   //
2456   // make a distortion map out ou fthe residual histogram
2457   // Results are written to the debug streamer - pcstream
2458   // Parameters:
2459   //   his0       - input (4D) residual histogram
2460   //   pcstream   - file to write the tree
2461   //   run        - run number
2462   //   refX       - track matching reference X
2463   //   type       - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
2464   // THnSparse axes:
2465   // OBJ: TAxis     #Delta  #Delta
2466   // OBJ: TAxis     tanTheta        tan(#Theta)
2467   // OBJ: TAxis     phi     #phi
2468   // OBJ: TAxis     snp     snp
2469
2470   // marian.ivanov@cern.ch
2471   const Int_t kMinEntries=10;
2472   Double_t bz=AliTrackerBase::GetBz();
2473   Int_t idim[4]={0,1,2,3};
2474   //
2475   //
2476   //
2477   Int_t nbins3=his0->GetAxis(3)->GetNbins();
2478   Int_t first3=his0->GetAxis(3)->GetFirst();
2479   Int_t last3 =his0->GetAxis(3)->GetLast();
2480   //
2481   for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){   // axis 3 - local angle
2482     his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3));
2483     Double_t      x3= his0->GetAxis(3)->GetBinCenter(ibin3);
2484     THnSparse * his3= his0->Projection(3,idim);         //projected histogram according selection 3
2485     //
2486     Int_t nbins2    = his3->GetAxis(2)->GetNbins();
2487     Int_t first2    = his3->GetAxis(2)->GetFirst();
2488     Int_t last2     = his3->GetAxis(2)->GetLast();
2489     //
2490     for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){   // axis 2 - phi
2491       his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2));
2492       Double_t      x2= his3->GetAxis(2)->GetBinCenter(ibin2);
2493       THnSparse * his2= his3->Projection(2,idim);         //projected histogram according selection 2
2494       Int_t nbins1     = his2->GetAxis(1)->GetNbins();
2495       Int_t first1     = his2->GetAxis(1)->GetFirst();
2496       Int_t last1      = his2->GetAxis(1)->GetLast();
2497       for (Int_t ibin1=first1; ibin1<last1; ibin1++){   //axis 1 - theta
2498         //
2499         Double_t       x1= his2->GetAxis(1)->GetBinCenter(ibin1);
2500         his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
2501         if (TMath::Abs(x1)<0.1){
2502           if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1));
2503           if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1));
2504         }
2505         if (TMath::Abs(x1)<0.06){
2506           his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
2507         }
2508         TH1 * hisDelta = his2->Projection(0);
2509         //
2510         Double_t entries = hisDelta->GetEntries();
2511         Double_t mean=0, rms=0;
2512         if (entries>kMinEntries){
2513           mean    = hisDelta->GetMean();
2514           rms = hisDelta->GetRMS();
2515         }
2516         Double_t sector = 9.*x2/TMath::Pi();
2517         if (sector<0) sector+=18;
2518         Double_t dsec = sector-Int_t(sector)-0.5;
2519         Double_t z=refX*x1;
2520         (*pcstream)<<hname<<
2521           "run="<<run<<
2522           "bz="<<bz<<
2523           "theta="<<x1<<
2524           "phi="<<x2<<
2525           "z="<<z<<            // dummy z
2526           "snp="<<x3<<
2527           "entries="<<entries<<
2528           "mean="<<mean<<
2529           "rms="<<rms<<
2530           "refX="<<refX<<   // track matching refernce plane
2531           "type="<<type<<   //
2532           "sector="<<sector<<
2533           "dsec="<<dsec<<
2534           "\n";
2535         delete hisDelta;
2536         //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
2537       }
2538       delete his2;
2539     }
2540     delete his3;
2541   }
2542 }
2543
2544
2545
2546
2547 void   AliTPCCorrection::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){
2548   //
2549   // make a distortion map out ou fthe residual histogram
2550   // Results are written to the debug streamer - pcstream
2551   // Parameters:
2552   //   his0       - input (4D) residual histogram
2553   //   pcstream   - file to write the tree
2554   //   run        - run number
2555   //   refX       - track matching reference X
2556   //   type       - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
2557   // marian.ivanov@cern.ch
2558   //
2559   //  Histo axeses
2560   //   Collection name='TObjArray', class='TObjArray', size=16
2561   //  0. OBJ: TAxis     #Delta  #Delta
2562   //  1. OBJ: TAxis     N_{cl}  N_{cl}
2563   //  2. OBJ: TAxis     dca_{r} (cm)    dca_{r} (cm)
2564   //  3. OBJ: TAxis     z (cm)  z (cm)
2565   //  4. OBJ: TAxis     sin(#phi)       sin(#phi)
2566   //  5. OBJ: TAxis     tan(#theta)     tan(#theta)
2567   //  6. OBJ: TAxis     1/pt (1/GeV)    1/pt (1/GeV)
2568   //  7. OBJ: TAxis     pt (GeV)        pt (GeV)
2569   //  8. OBJ: TAxis     alpha   alpha
2570   const Int_t kMinEntries=10;
2571   //
2572   //  1. make default selections
2573   //
2574   TH1 * hisDelta=0;
2575   Int_t idim0[4]={0 , 5, 8,  3};   // delta, theta, alpha, z
2576   hisInput->GetAxis(1)->SetRangeUser(110,190);   //long tracks
2577   hisInput->GetAxis(2)->SetRangeUser(-10,35);    //tracks close to beam pipe
2578   hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance
2579   hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks"
2580   hisDelta= hisInput->Projection(0);
2581   hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS());
2582   delete hisDelta;
2583   THnSparse *his0=  hisInput->Projection(4,idim0);
2584   //
2585   // 2. Get mean in diferent bins
2586   //
2587   Int_t nbins1=his0->GetAxis(1)->GetNbins();
2588   Int_t first1=his0->GetAxis(1)->GetFirst();
2589   Int_t last1 =his0->GetAxis(1)->GetLast();
2590   //
2591   Double_t bz=AliTrackerBase::GetBz();
2592   Int_t idim[4]={0,1, 2,  3};  // delta, theta,alpha,z
2593   //
2594   for (Int_t ibin1=first1; ibin1<=last1; ibin1++){   //axis 1 - theta
2595     //
2596     Double_t       x1= his0->GetAxis(1)->GetBinCenter(ibin1);
2597     his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
2598     //
2599     THnSparse * his1 = his0->Projection(4,idim);  // projected histogram according range1
2600     Int_t nbins3     = his1->GetAxis(3)->GetNbins();
2601     Int_t first3     = his1->GetAxis(3)->GetFirst();
2602     Int_t last3      = his1->GetAxis(3)->GetLast();
2603     //
2604     for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){   // axis 3 - z at "vertex"
2605       his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
2606       Double_t      x3= his1->GetAxis(3)->GetBinCenter(ibin3);
2607       if (ibin3<first3) {
2608         his1->GetAxis(3)->SetRangeUser(-1,1);
2609         x3=0;
2610       }
2611       THnSparse * his3= his1->Projection(4,idim);         //projected histogram according selection 3
2612       Int_t nbins2    = his3->GetAxis(2)->GetNbins();
2613       Int_t first2    = his3->GetAxis(2)->GetFirst();
2614       Int_t last2     = his3->GetAxis(2)->GetLast();
2615       //
2616       for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){
2617         his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
2618         Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
2619         hisDelta = his3->Projection(0);
2620         //
2621         Double_t entries = hisDelta->GetEntries();
2622         Double_t mean=0, rms=0;
2623         if (entries>kMinEntries){
2624           mean    = hisDelta->GetMean();
2625           rms = hisDelta->GetRMS();
2626         }
2627         Double_t sector = 9.*x2/TMath::Pi();
2628         if (sector<0) sector+=18;
2629         Double_t dsec = sector-Int_t(sector)-0.5;
2630         Double_t snp=0;  // dummy snp - equal 0
2631         (*pcstream)<<hname<<
2632           "run="<<run<<
2633           "bz="<<bz<<            // magnetic field
2634           "theta="<<x1<<         // theta
2635           "phi="<<x2<<           // phi (alpha)
2636           "z="<<x3<<             // z at "vertex"
2637           "snp="<<snp<<          // dummy snp
2638           "entries="<<entries<<  // entries in bin
2639           "mean="<<mean<<        // mean
2640           "rms="<<rms<<
2641           "refX="<<refX<<        // track matching refernce plane
2642           "type="<<type<<        // parameter type
2643           "sector="<<sector<<    // sector
2644           "dsec="<<dsec<<        // dummy delta sector
2645           "\n";
2646         delete hisDelta;
2647         printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
2648       }
2649       delete his3;
2650     }
2651     delete his1;
2652   }
2653   delete his0;
2654 }
2655
2656
2657
2658 void   AliTPCCorrection::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type){
2659   //
2660   // make a distortion map out of the residual histogram
2661   // Results are written to the debug streamer - pcstream
2662   // Parameters:
2663   //   his0       - input (4D) residual histogram
2664   //   pcstream   - file to write the tree
2665   //   run        - run number
2666   //   type       - 0- y 1-z,2 -snp, 3-theta
2667   // marian.ivanov@cern.ch
2668
2669   //Collection name='TObjArray', class='TObjArray', size=16
2670   //0  OBJ: TAxis     delta   delta
2671   //1  OBJ: TAxis     phi     phi
2672   //2  OBJ: TAxis     localX  localX
2673   //3  OBJ: TAxis     kY      kY
2674   //4  OBJ: TAxis     kZ      kZ
2675   //5  OBJ: TAxis     is1     is1
2676   //6  OBJ: TAxis     is0     is0
2677   //7. OBJ: TAxis     z       z
2678   //8. OBJ: TAxis     IsPrimary       IsPrimary
2679
2680   const Int_t kMinEntries=10;
2681   THnSparse * hisSector0=0;
2682   TH1 * htemp=0;    // histogram to calculate mean value of parameter
2683   Double_t bz=AliTrackerBase::GetBz();
2684
2685   //
2686   // Loop over pair of sector:
2687   // isPrim         - 8  ==> 8
2688   // isec0          - 6  ==> 7
2689   //   isec1        - 5  ==> 6
2690   //     refX       - 2  ==> 5
2691   //
2692   //     phi        - 1  ==> 4
2693   //       z        - 7  ==> 3
2694   //         snp    - 3  ==> 2
2695   //           theta- 4  ==> 1
2696   //                  0  ==> 0;
2697   for (Int_t isec0=0; isec0<72; isec0++){
2698     Int_t index0[9]={0, 4, 3, 7, 1, 2, 5, 6,8}; //regroup indeces
2699     //
2700     //hisInput->GetAxis(8)->SetRangeUser(-0.1,0.4);  // select secondaries only ? - get out later ?
2701     hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1);
2702     hisSector0=hisInput->Projection(7,index0);
2703     //
2704     //
2705     for (Int_t isec1=isec0+1; isec1<72; isec1++){
2706       //if (isec1!=isec0+36) continue;
2707       if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue;
2708       printf("Sectors %d\t%d\n",isec1,isec0);
2709       hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1);
2710       TH1 * hisX=hisSector0->Projection(5);
2711       Double_t refX= hisX->GetMean();
2712       delete hisX;
2713       TH1 *hisDelta=hisSector0->Projection(0);
2714       Double_t dmean = hisDelta->GetMean();
2715       Double_t drms = hisDelta->GetRMS();
2716       hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms);
2717       delete hisDelta;
2718       //
2719       //  1. make default selections
2720       //
2721       Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi }
2722       THnSparse *hisSector1=  hisSector0->Projection(5,idim0);
2723       //
2724       // 2. Get mean in diferent bins
2725       //
2726       Int_t idim[5]={0, 1, 2,  3, 4};  // {delta, theta-1,snp-2 ,z-3, phi-4}
2727       //
2728       //      Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins();
2729       Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst();
2730       Int_t lastPhi =hisSector1->GetAxis(4)->GetLast();
2731       //
2732       for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=1){   //axis 4 - phi
2733         //
2734         // Phi loop
2735         //
2736         Double_t       xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi);
2737         Double_t psec    = (9*xPhi/TMath::Pi());
2738         if (psec<0) psec+=18;
2739         Bool_t isOK0=kFALSE;
2740         Bool_t isOK1=kFALSE;
2741         if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.)  isOK0=kTRUE;
2742         if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.)  isOK1=kTRUE;
2743         if (!isOK0) continue;
2744         if (!isOK1) continue;
2745         //
2746         hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi));
2747         if (isec1!=isec0+36) {
2748           hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi));
2749         }
2750         //
2751         htemp = hisSector1->Projection(4);
2752         xPhi=htemp->GetMean();
2753         delete htemp;
2754         THnSparse * hisPhi = hisSector1->Projection(4,idim);
2755         //Int_t nbinsZ     = hisPhi->GetAxis(3)->GetNbins();
2756         Int_t firstZ     = hisPhi->GetAxis(3)->GetFirst();
2757         Int_t lastZ      = hisPhi->GetAxis(3)->GetLast();
2758         //
2759         for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=1){   // axis 3 - z
2760           //
2761           // Z loop
2762           //
2763           hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ));
2764           if (isec1!=isec0+36) {
2765             hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ));
2766           }
2767           htemp = hisPhi->Projection(3);
2768           Double_t      xZ= htemp->GetMean();
2769           delete htemp;
2770           THnSparse * hisZ= hisPhi->Projection(3,idim);
2771           //projected histogram according selection 3 -z
2772           //
2773           //
2774           //Int_t nbinsSnp    = hisZ->GetAxis(2)->GetNbins();
2775           Int_t firstSnp    = hisZ->GetAxis(2)->GetFirst();
2776           Int_t lastSnp     = hisZ->GetAxis(2)->GetLast();
2777           for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){   // axis 2 - snp
2778             //
2779             // Snp loop
2780             //
2781             hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp));
2782             if (isec1!=isec0+36) {
2783               hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp));
2784             }
2785             htemp = hisZ->Projection(2);
2786             Double_t      xSnp= htemp->GetMean();
2787             delete htemp;
2788             THnSparse * hisSnp= hisZ->Projection(2,idim);
2789             //projected histogram according selection 2 - snp
2790
2791             //Int_t nbinsTheta    = hisSnp->GetAxis(1)->GetNbins();
2792             Int_t firstTheta    = hisSnp->GetAxis(1)->GetFirst();
2793             Int_t lastTheta     = hisSnp->GetAxis(1)->GetLast();
2794             //
2795             for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){  // axis1 theta
2796
2797
2798               hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta));
2799               if (isec1!=isec0+36) {
2800                  hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta));
2801               }
2802               htemp = hisSnp->Projection(1);
2803               Double_t xTheta=htemp->GetMean();
2804               delete htemp;
2805               hisDelta = hisSnp->Projection(0);
2806               //
2807               Double_t entries = hisDelta->GetEntries();
2808               Double_t mean=0, rms=0;
2809               if (entries>kMinEntries){
2810                 mean    = hisDelta->GetMean();
2811                 rms = hisDelta->GetRMS();
2812               }
2813               Double_t sector = 9.*xPhi/TMath::Pi();
2814               if (sector<0) sector+=18;
2815               Double_t dsec = sector-Int_t(sector)-0.5;
2816               Int_t dtype=1;  // TPC alignment type
2817               (*pcstream)<<hname<<
2818                 "run="<<run<<
2819                 "bz="<<bz<<             // magnetic field
2820                 "ptype="<<type<<         // parameter type
2821                 "dtype="<<dtype<<         // parameter type
2822                 "isec0="<<isec0<<       // sector 0
2823                 "isec1="<<isec1<<       // sector 1
2824                 "sector="<<sector<<     // sector as float
2825                 "dsec="<<dsec<<         // delta sector
2826                 //
2827                 "theta="<<xTheta<<      // theta
2828                 "phi="<<xPhi<<          // phi (alpha)
2829                 "z="<<xZ<<              // z
2830                 "snp="<<xSnp<<          // snp
2831                 //
2832                 "entries="<<entries<<  // entries in bin
2833                 "mean="<<mean<<        // mean
2834                 "rms="<<rms<<          // rms
2835                 "refX="<<refX<<        // track matching reference plane
2836                 "\n";
2837               delete hisDelta;
2838               printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\n",isec0, isec1, xPhi,xZ,xSnp, xTheta, entries,mean);
2839               //
2840             }//ibinTheta
2841             delete hisSnp;
2842           } //ibinSnp
2843           delete hisZ;
2844         }//ibinZ
2845         delete hisPhi;
2846       }//ibinPhi
2847       delete hisSector1;
2848     }//isec1
2849     delete hisSector0;
2850   }//isec0
2851 }
2852
2853
2854
2855
2856
2857
2858
2859 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
2860   //
2861   // Store object in the OCDB
2862   // By default the object is stored in the current directory
2863   // default comment consit of user name and the date
2864   //
2865   TString ocdbStorage="";
2866   ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
2867   AliCDBMetaData *metaData= new AliCDBMetaData();
2868   metaData->SetObjectClassName("AliTPCCorrection");
2869   metaData->SetResponsible("Marian Ivanov");
2870   metaData->SetBeamPeriod(1);
2871   metaData->SetAliRootVersion("05-25-01"); //root version
2872   TString userName=gSystem->GetFromPipe("echo $USER");
2873   TString date=gSystem->GetFromPipe("date");
2874
2875   if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
2876   if (comment) metaData->SetComment(comment);
2877   AliCDBId* id1=NULL;
2878   id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
2879   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
2880   gStorage->Put(this, (*id1), metaData);
2881 }
2882
2883
2884 void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
2885   //
2886   // Fast method to simulate the influence of the given distortion on the vertex reconstruction
2887   //
2888
2889   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
2890   if (!magF) AliError("Magneticd field - not initialized");
2891   Double_t bz = magF->SolenoidField(); //field in kGauss
2892   printf("bz: %f\n",bz);
2893   AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
2894
2895   TObjArray   aTrk;              // Original Track array of Aside
2896   TObjArray   daTrk;             // Distorted Track array of A side
2897   UShort_t    *aId = new UShort_t[nTracks];      // A side Track ID
2898   TObjArray   cTrk;
2899   TObjArray   dcTrk;
2900   UShort_t    *cId = new UShort_t [nTracks];
2901   Int_t id=0;
2902   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2903   TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
2904   fpt.SetParameters(7.24,0.120);
2905   fpt.SetNpx(10000);
2906   for(Int_t nt=0; nt<nTracks; nt++){
2907     Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
2908     Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
2909     Double_t pt = fpt.GetRandom(); // momentum for f1
2910     //   printf("phi %lf  eta %lf pt %lf\n",phi,eta,pt);
2911     Short_t sign=1;
2912     if(gRandom->Rndm() < 0.5){
2913       sign =1;
2914     }else{
2915       sign=-1;
2916     }
2917
2918     Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
2919     Double_t pxyz[3];
2920     pxyz[0]=pt*TMath::Cos(phi);
2921     pxyz[1]=pt*TMath::Sin(phi);
2922     pxyz[2]=pt*TMath::Tan(theta);
2923     Double_t cv[21]={0};
2924     AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
2925
2926     Double_t refX=1.;
2927     Int_t dir=-1;
2928     AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir,  NULL);
2929     if (!td) continue;
2930     if (pcstream) (*pcstream)<<"track"<<
2931       "eta="<<eta<<
2932       "theta="<<theta<<
2933       "tOrig.="<<t<<
2934       "td.="<<td<<
2935       "\n";
2936     if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
2937       if (td){
2938         daTrk.AddLast(td);
2939         aTrk.AddLast(t);
2940         Int_t nn=aTrk.GetEntriesFast();
2941         aId[nn]=id;
2942       }
2943     }else if(( eta<-0.07 )&&( eta>-etaCuts )){
2944       if (td){
2945         dcTrk.AddLast(td);
2946         cTrk.AddLast(t);
2947         Int_t nn=cTrk.GetEntriesFast();
2948         cId[nn]=id;
2949       }
2950     }
2951     id++;
2952   }// end of track loop
2953
2954   vertexer->SetTPCMode();
2955   vertexer->SetConstraintOff();
2956
2957   aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));
2958   avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
2959   cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));
2960   cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
2961   if (pcstream) (*pcstream)<<"vertex"<<
2962     "x="<<orgVertex[0]<<
2963     "y="<<orgVertex[1]<<
2964     "z="<<orgVertex[2]<<
2965     "av.="<<&aV<<              // distorted vertex A side
2966     "cv.="<<&cV<<              // distroted vertex C side
2967     "avO.="<<&avOrg<<         // original vertex A side
2968     "cvO.="<<&cvOrg<<
2969     "\n";
2970   delete []aId;
2971   delete []cId;
2972 }
2973
2974 void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
2975   //
2976   // make correction available for visualization using
2977   // TFormula, TFX and TTree::Draw
2978   // important in order to check corrections and also compute dervied variables
2979   // e.g correction partial derivatives
2980   //
2981   // NOTE - class is not owner of correction
2982   //
2983   if (!fgVisualCorrection) fgVisualCorrection=new TObjArray(10000);
2984   if (position>=fgVisualCorrection->GetEntriesFast())
2985     fgVisualCorrection->Expand((position+10)*2);
2986   fgVisualCorrection->AddAt(corr, position);
2987 }
2988
2989
2990
2991 Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t r, Double_t kZ, Int_t axisType, Int_t corrType){
2992   //
2993   // calculate the correction at given position - check the geffCorr
2994   //
2995   // corrType return values
2996   // 0 - delta R
2997   // 1 - delta RPhi
2998   // 2 - delta Z
2999   // 3 - delta RPHI
3000   //
3001   if (!fgVisualCorrection) return 0;
3002   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3003   if (!corr) return 0;
3004
3005   Double_t phi=sector*TMath::Pi()/9.;
3006   Double_t gx = r*TMath::Cos(phi);
3007   Double_t gy = r*TMath::Sin(phi);
3008   Double_t gz = r*kZ;
3009   Int_t nsector=(gz>=0) ? 0:18;
3010   //
3011   //
3012   //
3013   Float_t distPoint[3]={gx,gy,gz};
3014   corr->DistortPoint(distPoint, nsector);
3015   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3016   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3017   Double_t phi0=TMath::ATan2(gy,gx);
3018   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3019   if (axisType==0) return r1-r0;
3020   if (axisType==1) return (phi1-phi0)*r0;
3021   if (axisType==2) return distPoint[2]-gz;
3022   if (axisType==3) return (TMath::Cos(phi)*(distPoint[0]-gx)+ TMath::Cos(phi)*(distPoint[1]-gy));
3023   return phi1-phi0;
3024 }
3025
3026 Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
3027   //
3028   // return correction at given x,y,z
3029   //
3030   if (!fgVisualCorrection) return 0;
3031   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3032   if (!corr) return 0;
3033   Double_t phi0= TMath::ATan2(gy,gx);
3034   Int_t nsector=(gz>=0) ? 0:18;
3035   Float_t distPoint[3]={gx,gy,gz};
3036   corr->CorrectPoint(distPoint, nsector);
3037   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3038   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3039   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3040   if (axisType==0) return r1-r0;
3041   if (axisType==1) return (phi1-phi0)*r0;
3042   if (axisType==2) return distPoint[2]-gz;
3043   return phi1-phi0;
3044 }
3045
3046 Double_t AliTPCCorrection::GetCorrXYZDz(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType,Double_t delta){
3047   //
3048   // return correction at given x,y,z
3049   //
3050   if (!fgVisualCorrection) return 0;
3051   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3052   if (!corr) return 0;
3053   Double_t phi0= TMath::ATan2(gy,gx);
3054   Int_t nsector=(gz>=0) ? 0:18;
3055   Float_t distPoint[3]={gx,gy,gz};
3056   Float_t dxyz[3]={gx,gy,gz};
3057   //
3058   corr->GetCorrectionDz(distPoint, nsector,dxyz,delta);
3059   distPoint[0]+=dxyz[0];
3060   distPoint[1]+=dxyz[1];
3061   distPoint[2]+=dxyz[2];
3062   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3063   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3064   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3065   if (axisType==0) return r1-r0;
3066   if (axisType==1) return (phi1-phi0)*r0;
3067   if (axisType==2) return distPoint[2]-gz;
3068   return phi1-phi0;
3069 }
3070
3071 Double_t AliTPCCorrection::GetCorrXYZIntegrateZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType,Double_t delta){
3072   //
3073   // return correction at given x,y,z
3074   //
3075   if (!fgVisualCorrection) return 0;
3076   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3077   if (!corr) return 0;
3078   Double_t phi0= TMath::ATan2(gy,gx);
3079   Int_t nsector=(gz>=0) ? 0:18;
3080   Float_t distPoint[3]={gx,gy,gz};
3081   Float_t dxyz[3]={gx,gy,gz};
3082   //
3083   corr->GetCorrectionIntegralDz(distPoint, nsector,dxyz,delta);
3084   distPoint[0]+=dxyz[0];
3085   distPoint[1]+=dxyz[1];
3086   distPoint[2]+=dxyz[2];
3087   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3088   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3089   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3090   if (axisType==0) return r1-r0;
3091   if (axisType==1) return (phi1-phi0)*r0;
3092   if (axisType==2) return distPoint[2]-gz;
3093   return phi1-phi0;
3094 }
3095
3096
3097 Double_t AliTPCCorrection::GetDistXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
3098   //
3099   // return correction at given x,y,z
3100   //
3101   if (!fgVisualCorrection) return 0;
3102   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3103   if (!corr) return 0;
3104   Double_t phi0= TMath::ATan2(gy,gx);
3105   Int_t nsector=(gz>=0) ? 0:18;
3106   Float_t distPoint[3]={gx,gy,gz};
3107   corr->DistortPoint(distPoint, nsector);
3108   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3109   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3110   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3111   if (axisType==0) return r1-r0;
3112   if (axisType==1) return (phi1-phi0)*r0;
3113   if (axisType==2) return distPoint[2]-gz;
3114   return phi1-phi0;
3115 }
3116
3117 Double_t AliTPCCorrection::GetDistXYZDz(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType,Double_t delta){
3118   //
3119   // return correction at given x,y,z
3120   //
3121   if (!fgVisualCorrection) return 0;
3122   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3123   if (!corr) return 0;
3124   Double_t phi0= TMath::ATan2(gy,gx);
3125   Int_t nsector=(gz>=0) ? 0:18;
3126   Float_t distPoint[3]={gx,gy,gz};
3127   Float_t dxyz[3]={gx,gy,gz};
3128   //
3129   corr->GetDistortionDz(distPoint, nsector,dxyz,delta);
3130   distPoint[0]+=dxyz[0];
3131   distPoint[1]+=dxyz[1];
3132   distPoint[2]+=dxyz[2];
3133   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3134   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3135   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3136   if (axisType==0) return r1-r0;
3137   if (axisType==1) return (phi1-phi0)*r0;
3138   if (axisType==2) return distPoint[2]-gz;
3139   return phi1-phi0;
3140 }
3141
3142 Double_t AliTPCCorrection::GetDistXYZIntegrateZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType,Double_t delta){
3143   //
3144   // return correction at given x,y,z
3145   //
3146   if (!fgVisualCorrection) return 0;
3147   AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3148   if (!corr) return 0;
3149   Double_t phi0= TMath::ATan2(gy,gx);
3150   Int_t nsector=(gz>=0) ? 0:18;
3151   Float_t distPoint[3]={gx,gy,gz};
3152   Float_t dxyz[3]={gx,gy,gz};
3153   //
3154   corr->GetDistortionIntegralDz(distPoint, nsector,dxyz,delta);
3155   distPoint[0]+=dxyz[0];
3156   distPoint[1]+=dxyz[1];
3157   distPoint[2]+=dxyz[2];
3158   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3159   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3160   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3161   if (axisType==0) return r1-r0;
3162   if (axisType==1) return (phi1-phi0)*r0;
3163   if (axisType==2) return distPoint[2]-gz;
3164   return phi1-phi0;
3165 }
3166
3167
3168
3169 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){
3170   //
3171   // Make a laser fit tree for global minimization
3172   //
3173   AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
3174   AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
3175   if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
3176   correction->AddVisualCorrection(correction,0);  //register correction
3177
3178   //  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
3179   //AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
3180   //
3181   const Double_t cutErrY=0.05;
3182   const Double_t kSigmaCut=4;
3183   //  const Double_t cutErrZ=0.03;
3184   const Double_t kEpsilon=0.00000001;
3185   //  const Double_t kMaxDist=1.;  // max distance - space correction
3186   TVectorD *vecdY=0;
3187   TVectorD *vecdZ=0;
3188   TVectorD *veceY=0;
3189   TVectorD *veceZ=0;
3190   AliTPCLaserTrack *ltr=0;
3191   AliTPCLaserTrack::LoadTracks();
3192   tree->SetBranchAddress("dY.",&vecdY);
3193   tree->SetBranchAddress("dZ.",&vecdZ);
3194   tree->SetBranchAddress("eY.",&veceY);
3195   tree->SetBranchAddress("eZ.",&veceZ);
3196   tree->SetBranchAddress("LTr.",&ltr);
3197   Int_t entries= tree->GetEntries();
3198   TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
3199   Double_t bz=AliTrackerBase::GetBz();
3200   //
3201   //  Double_t globalXYZ[3];
3202   //Double_t globalXYZCorr[3];
3203   for (Int_t ientry=0; ientry<entries; ientry++){
3204     tree->GetEntry(ientry);
3205     if (!ltr->GetVecGX()){
3206       ltr->UpdatePoints();
3207     }
3208     //
3209     TVectorD fit10(5);
3210     TVectorD fit5(5);
3211     printf("Entry\t%d\n",ientry);
3212     for (Int_t irow0=0; irow0<158; irow0+=1){
3213       //
3214       TLinearFitter fitter10(4,"hyp3");
3215       TLinearFitter fitter5(2,"hyp1");
3216       Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0];
3217       if (sector<0) continue;
3218       //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue;
3219
3220       Double_t refX= (*ltr->GetVecLX())[irow0];
3221       Int_t firstRow1 = TMath::Max(irow0-10,0);
3222       Int_t lastRow1  = TMath::Min(irow0+10,158);
3223       Double_t padWidth=(irow0<64)?0.4:0.6;
3224       // make long range fit
3225       for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){
3226         if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
3227         if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
3228         if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
3229         Double_t idealX= (*ltr->GetVecLX())[irow1];
3230         Double_t idealY= (*ltr->GetVecLY())[irow1];
3231         //      Double_t idealZ= (*ltr->GetVecLZ())[irow1];
3232         Double_t gx= (*ltr->GetVecGX())[irow1];
3233         Double_t gy= (*ltr->GetVecGY())[irow1];
3234         Double_t gz= (*ltr->GetVecGZ())[irow1];
3235         Double_t measY=(*vecdY)[irow1]+idealY;
3236         Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
3237         // deltaR = R distorted -R ideal
3238         Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
3239         fitter10.AddPoint(xxx,measY,1);
3240       }
3241       Bool_t isOK=kTRUE;
3242       Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
3243       Double_t mean10  =0;//   fitter10.GetParameter(0);
3244       Double_t slope10  =0;//   fitter10.GetParameter(0);
3245       Double_t cosPart10  = 0;//  fitter10.GetParameter(2);
3246       Double_t sinPart10   =0;//  fitter10.GetParameter(3);
3247
3248       if (fitter10.GetNpoints()>10){
3249         fitter10.Eval();
3250         rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
3251         mean10      =   fitter10.GetParameter(0);
3252         slope10     =   fitter10.GetParameter(1);
3253         cosPart10   =   fitter10.GetParameter(2);
3254         sinPart10   =  fitter10.GetParameter(3);
3255         //
3256         // make short range fit
3257         //
3258         for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){
3259           if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
3260           if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
3261           if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
3262           Double_t idealX= (*ltr->GetVecLX())[irow1];
3263           Double_t idealY= (*ltr->GetVecLY())[irow1];
3264           //      Double_t idealZ= (*ltr->GetVecLZ())[irow1];
3265           Double_t gx= (*ltr->GetVecGX())[irow1];
3266           Double_t gy= (*ltr->GetVecGY())[irow1];
3267           Double_t gz= (*ltr->GetVecGZ())[irow1];
3268           Double_t measY=(*vecdY)[irow1]+idealY;
3269           Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
3270           // deltaR = R distorted -R ideal
3271           Double_t expY= mean10+slope10*(idealX+deltaR-refX);
3272           if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue;
3273           //
3274           Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
3275           Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
3276           fitter5.AddPoint(xxx,measY-corr,1);
3277         }
3278       }else{
3279         isOK=kFALSE;
3280       }
3281       if (fitter5.GetNpoints()<8) isOK=kFALSE;
3282
3283       Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
3284       Double_t offset5  =0;//  fitter5.GetParameter(0);
3285       Double_t slope5   =0;//  fitter5.GetParameter(0);
3286       if (isOK){
3287         fitter5.Eval();
3288         rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
3289         offset5  =  fitter5.GetParameter(0);
3290         slope5   =  fitter5.GetParameter(0);
3291       }
3292       //
3293       Double_t dtype=5;
3294       Double_t ptype=0;
3295       Double_t phi   =(*ltr->GetVecPhi())[irow0];
3296       Double_t theta =ltr->GetTgl();
3297       Double_t mean=(vecdY)->GetMatrixArray()[irow0];
3298       Double_t gx=0,gy=0,gz=0;
3299       Double_t snp = (*ltr->GetVecP2())[irow0];
3300       Int_t bundle= ltr->GetBundle();
3301       Int_t id= ltr->GetId();
3302       //      Double_t rms = err->GetMatrixArray()[irow];
3303       //
3304       gx = (*ltr->GetVecGX())[irow0];
3305       gy = (*ltr->GetVecGY())[irow0];
3306       gz = (*ltr->GetVecGZ())[irow0];
3307       Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0);
3308       fitter10.GetParameters(fit10);
3309       fitter5.GetParameters(fit5);
3310       Double_t idealY= (*ltr->GetVecLY())[irow0];
3311       Double_t measY=(*vecdY)[irow0]+idealY;
3312       Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
3313       if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE;
3314       //
3315       (*pcstream)<<"fitFull"<<  // dumpe also intermediate results
3316         "bz="<<bz<<         // magnetic filed used
3317         "dtype="<<dtype<<   // detector match type
3318         "ptype="<<ptype<<   // parameter type
3319         "theta="<<theta<<   // theta
3320         "phi="<<phi<<       // phi
3321         "snp="<<snp<<       // snp
3322         "sector="<<sector<<
3323         "bundle="<<bundle<<
3324 //      //      "dsec="<<dsec<<
3325         "refX="<<refX<<      // reference radius
3326         "gx="<<gx<<         // global position
3327         "gy="<<gy<<         // global position
3328         "gz="<<gz<<         // global position
3329         "dRrec="<<dRrec<<      // delta Radius in reconstruction
3330         "id="<<id<<     //bundle
3331         "rms10="<<rms10<<
3332         "rms5="<<rms5<<
3333         "fit10.="<<&fit10<<
3334         "fit5.="<<&fit5<<
3335         "measY="<<measY<<
3336         "mean="<<mean<<
3337         "idealY="<<idealY<<
3338         "corr="<<corr<<
3339         "isOK="<<isOK<<
3340         "\n";
3341     }
3342   }
3343   delete pcstream;
3344 }