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