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