TPC/AliTPCCorrection.cxx

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