adding the tansformation names to the trees
[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
42 #include <TH2F.h>
43 #include <TMath.h>
44 #include <TROOT.h>
45 #include <TTreeStream.h>
46
47 #include  "AliExternalTrackParam.h"
48 #include  "AliTrackPointArray.h"
49 #include  "TDatabasePDG.h"
50 #include  "AliTrackerBase.h"
51 #include  "AliTPCROC.h"
52
53
54 #include "AliTPCCorrection.h"
55
56 ClassImp(AliTPCCorrection)
57
58 // FIXME: the following values should come from the database
59 const Double_t AliTPCCorrection::fgkTPC_Z0   =249.7;     // nominal gating grid position 
60 const Double_t AliTPCCorrection::fgkIFCRadius= 83.06;    // Mean Radius of the Inner Field Cage ( 82.43 min,  83.70 max) (cm)
61 const Double_t AliTPCCorrection::fgkOFCRadius=254.5;     // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
62 const Double_t AliTPCCorrection::fgkZOffSet  = 0.2;      // Offset from CE: calculate all distortions closer to CE as if at this point
63 const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
64 const Double_t AliTPCCorrection::fgkGG       =-70.0;     // Gating Grid voltage (volts)
65
66
67 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
68 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] =  {   
69 84.0,   84.5,   85.0,   85.5,   86.0,  87.0,    88.0,
70 90.0,   92.0,   94.0,   96.0,   98.0,  100.0,  102.0,  104.0,  106.0,  108.0, 
71 110.0,  112.0,  114.0,  116.0,  118.0,  120.0,  122.0,  124.0,  126.0,  128.0, 
72 130.0,  132.0,  134.0,  136.0,  138.0,  140.0,  142.0,  144.0,  146.0,  148.0, 
73 150.0,  152.0,  154.0,  156.0,  158.0,  160.0,  162.0,  164.0,  166.0,  168.0, 
74 170.0,  172.0,  174.0,  176.0,  178.0,  180.0,  182.0,  184.0,  186.0,  188.0,
75 190.0,  192.0,  194.0,  196.0,  198.0,  200.0,  202.0,  204.0,  206.0,  208.0,
76 210.0,  212.0,  214.0,  216.0,  218.0,  220.0,  222.0,  224.0,  226.0,  228.0,
77 230.0,  232.0,  234.0,  236.0,  238.0,  240.0,  242.0,  244.0,  246.0,  248.0,
78 249.0,  249.5,  250.0,  251.5,  252.0  } ;
79   
80 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ]     =   { 
81 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
82 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
83 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
84 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
85 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
86 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
87 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
88 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
89 -100.0,  -98.0,  -96.0,  -94.0,  -92.0,  -90.0,  -88.0,  -86.0,  -84.0,  -82.0,
90 -80.0,  -78.0,  -76.0,  -74.0,  -72.0,  -70.0,  -68.0,  -66.0,  -64.0,  -62.0,
91 -60.0,  -58.0,  -56.0,  -54.0,  -52.0,  -50.0,  -48.0,  -46.0,  -44.0,  -42.0,
92 -40.0,  -38.0,  -36.0,  -34.0,  -32.0,  -30.0,  -28.0,  -26.0,  -24.0,  -22.0,
93 -20.0,  -18.0,  -16.0,  -14.0,  -12.0,  -10.0,   -8.0,   -6.0,   -4.0,   -2.0,
94 -1.0,   -0.5,   -0.2,   -0.1,  -0.05,   0.05,    0.1,    0.2,    0.5,    1.0, 
95  2.0,    4.0,    6.0,    8.0,   10.0,   12.0,   14.0,   16.0,   18.0,   20.0,
96  22.0,   24.0,   26.0,   28.0,   30.0,   32.0,   34.0,   36.0,   38.0,   40.0, 
97  42.0,   44.0,   46.0,   48.0,   50.0,   52.0,   54.0,   56.0,   58.0,   60.0, 
98  62.0,   64.0,   66.0,   68.0,   70.0,   72.0,   74.0,   76.0,   78.0,   80.0, 
99  82.0,   84.0,   86.0,   88.0,   90.0,   92.0,   94.0,   96.0,   98.0,  100.0, 
100 102.0,  104.0,  106.0,  108.0,  110.0,  112.0,  114.0,  116.0,  118.0,  120.0, 
101 122.0,  124.0,  126.0,  128.0,  130.0,  132.0,  134.0,  136.0,  138.0,  140.0, 
102 142.0,  144.0,  146.0,  148.0,  150.0,  152.0,  154.0,  156.0,  158.0,  160.0, 
103 162.0,  164.0,  166.0,  168.0,  170.0,  172.0,  174.0,  176.0,  178.0,  180.0, 
104 182.0,  184.0,  186.0,  188.0,  190.0,  192.0,  194.0,  196.0,  198.0,  200.0,
105 202.0,  204.0,  206.0,  208.0,  210.0,  212.0,  214.0,  216.0,  218.0,  220.0,
106 222.0,  224.0,  226.0,  228.0,  230.0,  232.0,  234.0,  236.0,  238.0,  240.0,
107 242.0,  243.0,  244.0,  245.0,  246.0,  247.0,  248.0,  248.5,  249.0,  249.5   } ;
108
109
110
111 AliTPCCorrection::AliTPCCorrection() 
112   : TNamed("correction_unity","unity"),fJLow(0),fKLow(0)
113 {
114   //
115   // default constructor
116   //
117 }
118
119 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
120   : TNamed(name,title),fJLow(0),fKLow(0)
121 {
122   //
123   // default constructor, that set the name and title
124   //
125 }
126
127 AliTPCCorrection::~AliTPCCorrection() {
128   // 
129   // virtual destructor
130   //
131 }
132
133 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
134   //
135   // Corrects the initial coordinates x (cartesian coordinates)
136   // according to the given effect (inherited classes)
137   // roc represents the TPC read out chamber (offline numbering convention)
138   //
139   Float_t dx[3];
140   GetCorrection(x,roc,dx);
141   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
142 }
143
144 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
145   //
146   // Corrects the initial coordinates x (cartesian coordinates) and stores the new 
147   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
148   // roc represents the TPC read out chamber (offline numbering convention)
149   //
150   Float_t dx[3];
151   GetCorrection(x,roc,dx);
152   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
153 }
154
155 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
156   //
157   // Distorts the initial coordinates x (cartesian coordinates)
158   // according to the given effect (inherited classes)
159   // roc represents the TPC read out chamber (offline numbering convention)
160   //
161   Float_t dx[3];
162   GetDistortion(x,roc,dx);
163   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
164 }
165
166 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
167   //
168   // Distorts the initial coordinates x (cartesian coordinates) and stores the new 
169   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
170   // roc represents the TPC read out chamber (offline numbering convention)
171   //
172   Float_t dx[3];
173   GetDistortion(x,roc,dx);
174   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
175 }
176
177 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
178   //
179   // This function delivers the correction values dx in respect to the inital coordinates x
180   // roc represents the TPC read out chamber (offline numbering convention)
181   // Note: The dx is overwritten by the inherited effectice class ...
182   //
183   for (Int_t j=0;j<3;++j) { dx[j]=0.; }
184 }
185
186 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
187   //
188   // This function delivers the distortion values dx in respect to the inital coordinates x
189   // roc represents the TPC read out chamber (offline numbering convention)
190   //
191   GetCorrection(x,roc,dx);
192   for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
193 }
194
195 void AliTPCCorrection::Init() {
196   //
197   // Initialization funtion (not used at the moment)
198   //
199 }
200
201 void AliTPCCorrection::Print(Option_t* /*option*/) const {
202   //
203   // Print function to check which correction classes are used 
204   // option=="d" prints details regarding the setted magnitude 
205   // option=="a" prints the C0 and C1 coefficents for calibration purposes
206   //
207   printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
208 }
209
210 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t /*t1*/,Float_t /*t2*/) {
211   //
212   // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
213   // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
214   // calibration run
215   //
216   // SetOmegaTauT1T2(omegaTau, t1, t2);
217 }
218
219 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
220   //
221   // Simple plot functionality.
222   // Returns a 2d hisogram which represents the corrections in radial direction (dr)
223   // in respect to position z within the XY plane.
224   // The histogramm has nx times ny entries. 
225   //
226   
227   TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
228                      nx,-250.,250.,ny,-250.,250.);
229   Float_t x[3],dx[3];
230   x[2]=z;
231   Int_t roc=z>0.?0:18; // FIXME
232   for (Int_t iy=1;iy<=ny;++iy) {
233     x[1]=h->GetYaxis()->GetBinCenter(iy);
234     for (Int_t ix=1;ix<=nx;++ix) {
235       x[0]=h->GetXaxis()->GetBinCenter(ix);
236       GetCorrection(x,roc,dx);
237       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
238       if (90.<=r0 && r0<=250.) {
239         Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
240         h->SetBinContent(ix,iy,r1-r0);
241       }
242       else
243         h->SetBinContent(ix,iy,0.);
244     }
245   }
246   return h;
247 }
248
249 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
250   //
251   // Simple plot functionality.
252   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi) 
253   // in respect to position z within the XY plane.
254   // The histogramm has nx times ny entries. 
255   //
256
257   TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
258                      nx,-250.,250.,ny,-250.,250.);
259   Float_t x[3],dx[3];
260   x[2]=z;
261   Int_t roc=z>0.?0:18; // FIXME
262   for (Int_t iy=1;iy<=ny;++iy) {
263     x[1]=h->GetYaxis()->GetBinCenter(iy);
264     for (Int_t ix=1;ix<=nx;++ix) {
265       x[0]=h->GetXaxis()->GetBinCenter(ix);
266       GetCorrection(x,roc,dx);
267       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
268       if (90.<=r0 && r0<=250.) {
269         Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
270         Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
271
272         Float_t dphi=phi1-phi0;
273         if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
274         if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
275       
276         h->SetBinContent(ix,iy,r0*dphi);
277       }
278       else
279         h->SetBinContent(ix,iy,0.);
280     }
281   }
282   return h;
283 }
284
285 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
286   //
287   // Simple plot functionality.
288   // Returns a 2d hisogram which represents the corrections in r direction (dr) 
289   // in respect to angle phi within the ZR plane.
290   // The histogramm has nx times ny entries. 
291   //
292   TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
293                      nz,-250.,250.,nr,85.,250.);
294   Float_t x[3],dx[3];
295   for (Int_t ir=1;ir<=nr;++ir) {
296     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
297     x[0]=radius*TMath::Cos(phi);
298     x[1]=radius*TMath::Sin(phi);
299     for (Int_t iz=1;iz<=nz;++iz) {
300       x[2]=h->GetXaxis()->GetBinCenter(iz);
301       Int_t roc=x[2]>0.?0:18; // FIXME
302       GetCorrection(x,roc,dx);
303       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
304       Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
305       h->SetBinContent(iz,ir,r1-r0);
306     }
307   }
308   printf("SDF\n");
309   return h;
310
311 }
312
313 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
314   //
315   // Simple plot functionality.
316   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi) 
317   // in respect to angle phi within the ZR plane.
318   // The histogramm has nx times ny entries. 
319   //
320   TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
321                      nz,-250.,250.,nr,85.,250.);
322   Float_t x[3],dx[3];
323   for (Int_t iz=1;iz<=nz;++iz) {
324     x[2]=h->GetXaxis()->GetBinCenter(iz);
325     Int_t roc=x[2]>0.?0:18; // FIXME
326     for (Int_t ir=1;ir<=nr;++ir) {
327       Float_t radius=h->GetYaxis()->GetBinCenter(ir);
328       x[0]=radius*TMath::Cos(phi);
329       x[1]=radius*TMath::Sin(phi);
330       GetCorrection(x,roc,dx);
331       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
332       Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
333       Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
334       
335       Float_t dphi=phi1-phi0;
336       if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
337       if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
338       
339       h->SetBinContent(iz,ir,r0*dphi);
340     }
341   }
342   return h;
343 }
344
345 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
346                                    const char *xlabel,const char *ylabel,const char *zlabel,
347                                   Int_t nbinsx,Double_t xlow,Double_t xup,
348                                   Int_t nbinsy,Double_t ylow,Double_t yup) {
349   //
350   // Helper function to create a 2d histogramm of given size
351   //
352   
353   TString hname=name;
354   Int_t i=0;
355   if (gDirectory) {
356     while (gDirectory->FindObject(hname.Data())) {
357       hname =name;
358       hname+="_";
359       hname+=i;
360       ++i;
361     }
362   }
363   TH2F *h=new TH2F(hname.Data(),title,
364                    nbinsx,xlow,xup,
365                    nbinsy,ylow,yup);
366   h->GetXaxis()->SetTitle(xlabel);
367   h->GetYaxis()->SetTitle(ylabel);
368   h->GetZaxis()->SetTitle(zlabel);
369   h->SetStats(0);
370   return h;
371 }
372
373
374 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
375
376 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z, 
377                                                   const Double_t er[kNZ][kNR], Double_t &er_value )
378 {
379   //
380   // Interpolate table - 2D interpolation
381   //
382   Double_t save_er[10] ;
383
384   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
385   Search( kNR,   fgkRList,  r,   fKLow   ) ;
386   if ( fJLow < 0 ) fJLow = 0 ;   // check if out of range
387   if ( fKLow < 0 ) fKLow = 0 ;
388   if ( fJLow + order  >=    kNZ - 1 ) fJLow =   kNZ - 1 - order ;
389   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
390
391   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
392       save_er[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
393   }
394   er_value = Interpolate( &fgkZList[fJLow], save_er, order, z )   ;
395
396 }
397
398
399 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[], 
400                                        const Int_t order, const Double_t x )
401 {
402   //
403   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
404   //
405
406   Double_t y ;
407   if ( order == 2 ) {                // Quadratic Interpolation = 2 
408     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ; 
409     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ; 
410     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ; 
411   } else {                           // Linear Interpolation = 1
412     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
413   }
414
415   return (y);
416
417 }
418
419
420 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low )
421 {
422   //
423   // Search an ordered table by starting at the most recently used point
424   //
425
426   Long_t middle, high ;
427   Int_t  ascend = 0, increment = 1 ;
428
429   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
430   
431   if ( low < 0 || low > n-1 ) { 
432     low = -1 ; high = n ; 
433   } else {                                            // Ordered Search phase
434     if ( (Int_t)( x >= xArray[low] ) == ascend )  {
435       if ( low == n-1 ) return ;          
436       high = low + 1 ;
437       while ( (Int_t)( x >= xArray[high] ) == ascend ) {
438         low = high ;
439         increment *= 2 ;
440         high = low + increment ;
441         if ( high > n-1 )  {  high = n ; break ;  }
442       }
443     } else {
444       if ( low == 0 )  {  low = -1 ;  return ;  }
445       high = low - 1 ;
446       while ( (Int_t)( x < xArray[low] ) == ascend ) {
447         high = low ;
448         increment *= 2 ;
449         if ( increment >= high )  {  low = -1 ;  break ;  }
450         else  low = high - increment ;
451       }
452     }
453   }
454   
455   while ( (high-low) != 1 ) {                     // Binary Search Phase
456     middle = ( high + low ) / 2 ;
457     if ( (Int_t)( x >= xArray[middle] ) == ascend )
458       low = middle ;
459     else
460       high = middle ;
461   }
462   
463   if ( x == xArray[n-1] ) low = n-2 ;
464   if ( x == xArray[0]   ) low = 0 ;
465   
466 }
467
468
469 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(const AliExternalTrackParam * trackIn, Double_t refX, Int_t dir,TTreeSRedirector *pcstream){
470   //
471   // Fit the track parameters - without and with distortion
472   // 1. Space points in the TPC are simulated along the trajectory  
473   // 2. Space points distorted
474   // 3. Fits  the non distorted and distroted track to the reference plane at refX
475   // 4. For visualization and debugging  purposes the space points and tracks can be stored  in the tree - using the TTreeSRedirector functionality 
476   //
477   // trackIn   - input track parameters
478   // refX     - reference X to fit the track
479   // dir      - direction - out=1 or in=-1
480   // pcstream -  debug streamer to check the results
481   //
482   // see AliExternalTrackParam.h documentation:
483   // track1.fP[0] - local y (rphi)
484   // track1.fP[1] - z 
485   // track1.fP[2] - sinus of local inclination angle
486   // track1.fP[3] - tangent of deep angle
487   // track1.fP[4] - 1/pt
488   AliTPCROC * roc = AliTPCROC::Instance();
489   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
490   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
491   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
492     
493   const Double_t kMaxSnp = 0.85;  
494   const Double_t kSigmaY=0.1;
495   const Double_t kSigmaZ=0.1;
496   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
497
498   AliExternalTrackParam  track(*trackIn); // 
499   // generate points
500   AliTrackPointArray pointArray0(npoints0);
501   AliTrackPointArray pointArray1(npoints0);
502   Double_t xyz[3];
503   AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
504   //
505   // simulate the track
506   Int_t npoints=0;
507   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
508   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
509     AliTrackerBase::PropagateTrackTo(&track,radius,kMass,3,kTRUE,kMaxSnp);
510     track.GetXYZ(xyz);
511     AliTrackPoint pIn0;                               // space point          
512     AliTrackPoint pIn1;
513     Int_t roc= (xyz[2]>0)? 0:18;
514     pointArray0.GetPoint(pIn0,npoints);
515     pointArray1.GetPoint(pIn1,npoints);
516     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
517     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
518     DistortPoint(distPoint, roc);
519     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
520     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
521     //
522     track.Rotate(alpha);
523     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
524     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
525     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
526     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
527     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
528     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
529     pointArray0.AddPoint(npoints, &pIn0);      
530     pointArray1.AddPoint(npoints, &pIn1);
531     npoints++;
532     if (npoints>=npoints0) break;
533   }
534   //
535   // refit track
536   //
537   AliExternalTrackParam *track0=0;
538   AliExternalTrackParam *track1=0;
539   AliTrackPoint   point1,point2,point3;
540   if (dir==1) {  //make seed inner
541     pointArray0.GetPoint(point1,1);
542     pointArray0.GetPoint(point2,10);
543     pointArray0.GetPoint(point3,20);
544   }
545   if (dir==-1){ //make seed outer
546     pointArray0.GetPoint(point1,npoints-20);
547     pointArray0.GetPoint(point2,npoints-10);
548     pointArray0.GetPoint(point3,npoints-1);
549   }  
550   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
551   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
552
553
554   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
555     Int_t ipoint= (dir>0) ? ipoint: npoints-1-jpoint;
556     //
557     AliTrackPoint pIn0;
558     AliTrackPoint pIn1;
559     pointArray0.GetPoint(pIn0,ipoint);
560     pointArray1.GetPoint(pIn1,ipoint);
561     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
562     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
563     //
564     AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,1,kFALSE,kMaxSnp);
565     AliTrackerBase::PropagateTrackTo(track1,prot1.GetX(),kMass,1,kFALSE,kMaxSnp);
566     track.GetXYZ(xyz);
567     //
568     Double_t pointPos[2]={0,0};
569     Double_t pointCov[3]={0,0,0};
570     pointPos[0]=prot0.GetY();//local y
571     pointPos[1]=prot0.GetZ();//local z
572     pointCov[0]=prot0.GetCov()[3];//simay^2
573     pointCov[1]=prot0.GetCov()[4];//sigmayz
574     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
575     track0->Update(pointPos,pointCov);
576     //
577     pointPos[0]=prot1.GetY();//local y
578     pointPos[1]=prot1.GetZ();//local z
579     pointCov[0]=prot1.GetCov()[3];//simay^2
580     pointCov[1]=prot1.GetCov()[4];//sigmayz
581     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
582     track1->Update(pointPos,pointCov);
583   }
584
585   AliTrackerBase::PropagateTrackTo(track0,refX,kMass,2.,kTRUE,kMaxSnp);
586   track1->Rotate(track0->GetAlpha());
587   AliTrackerBase::PropagateTrackTo(track1,refX,kMass,2.,kFALSE,kMaxSnp);
588
589   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
590     "point0.="<<&pointArray0<<   //  points
591     "point1.="<<&pointArray1<<   //  distorted points
592     "trackIn.="<<&track<<       //  original track
593     "track0.="<<track0<<         //  fitted track
594     "track1.="<<track1<<         //  fitted distorted track
595     "\n";
596   delete track0;
597   return track1;
598 }
599
600