Fixing coding violation
[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
54 #include  "TRandom.h"
55 #include  "AliExternalTrackParam.h"
56 #include  "AliTrackPointArray.h"
57 #include  "TDatabasePDG.h"
58 #include  "AliTrackerBase.h"
59 #include  "AliTPCROC.h"
60 #include  "THnSparse.h"
61
62
63 #include "AliTPCCorrection.h"
64
65 ClassImp(AliTPCCorrection)
66
67 // FIXME: the following values should come from the database
68 const Double_t AliTPCCorrection::fgkTPCZ0    =249.7;     // nominal gating grid position 
69 const Double_t AliTPCCorrection::fgkIFCRadius= 83.06;    // Mean Radius of the Inner Field Cage ( 82.43 min,  83.70 max) (cm)
70 const Double_t AliTPCCorrection::fgkOFCRadius=254.5;     // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
71 const Double_t AliTPCCorrection::fgkZOffSet  = 0.2;      // Offset from CE: calculate all distortions closer to CE as if at this point
72 const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
73 const Double_t AliTPCCorrection::fgkGG       =-70.0;     // Gating Grid voltage (volts)
74
75
76 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
77 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] =  {   
78 84.0,   84.5,   85.0,   85.5,   86.0,  87.0,    88.0,
79 90.0,   92.0,   94.0,   96.0,   98.0,  100.0,  102.0,  104.0,  106.0,  108.0, 
80 110.0,  112.0,  114.0,  116.0,  118.0,  120.0,  122.0,  124.0,  126.0,  128.0, 
81 130.0,  132.0,  134.0,  136.0,  138.0,  140.0,  142.0,  144.0,  146.0,  148.0, 
82 150.0,  152.0,  154.0,  156.0,  158.0,  160.0,  162.0,  164.0,  166.0,  168.0, 
83 170.0,  172.0,  174.0,  176.0,  178.0,  180.0,  182.0,  184.0,  186.0,  188.0,
84 190.0,  192.0,  194.0,  196.0,  198.0,  200.0,  202.0,  204.0,  206.0,  208.0,
85 210.0,  212.0,  214.0,  216.0,  218.0,  220.0,  222.0,  224.0,  226.0,  228.0,
86 230.0,  232.0,  234.0,  236.0,  238.0,  240.0,  242.0,  244.0,  246.0,  248.0,
87 249.0,  249.5,  250.0,  251.5,  252.0  } ;
88   
89 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ]     =   { 
90 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
91 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
92 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
93 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
94 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
95 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
96 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
97 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
98 -100.0,  -98.0,  -96.0,  -94.0,  -92.0,  -90.0,  -88.0,  -86.0,  -84.0,  -82.0,
99 -80.0,  -78.0,  -76.0,  -74.0,  -72.0,  -70.0,  -68.0,  -66.0,  -64.0,  -62.0,
100 -60.0,  -58.0,  -56.0,  -54.0,  -52.0,  -50.0,  -48.0,  -46.0,  -44.0,  -42.0,
101 -40.0,  -38.0,  -36.0,  -34.0,  -32.0,  -30.0,  -28.0,  -26.0,  -24.0,  -22.0,
102 -20.0,  -18.0,  -16.0,  -14.0,  -12.0,  -10.0,   -8.0,   -6.0,   -4.0,   -2.0,
103 -1.0,   -0.5,   -0.2,   -0.1,  -0.05,   0.05,    0.1,    0.2,    0.5,    1.0, 
104  2.0,    4.0,    6.0,    8.0,   10.0,   12.0,   14.0,   16.0,   18.0,   20.0,
105  22.0,   24.0,   26.0,   28.0,   30.0,   32.0,   34.0,   36.0,   38.0,   40.0, 
106  42.0,   44.0,   46.0,   48.0,   50.0,   52.0,   54.0,   56.0,   58.0,   60.0, 
107  62.0,   64.0,   66.0,   68.0,   70.0,   72.0,   74.0,   76.0,   78.0,   80.0, 
108  82.0,   84.0,   86.0,   88.0,   90.0,   92.0,   94.0,   96.0,   98.0,  100.0, 
109 102.0,  104.0,  106.0,  108.0,  110.0,  112.0,  114.0,  116.0,  118.0,  120.0, 
110 122.0,  124.0,  126.0,  128.0,  130.0,  132.0,  134.0,  136.0,  138.0,  140.0, 
111 142.0,  144.0,  146.0,  148.0,  150.0,  152.0,  154.0,  156.0,  158.0,  160.0, 
112 162.0,  164.0,  166.0,  168.0,  170.0,  172.0,  174.0,  176.0,  178.0,  180.0, 
113 182.0,  184.0,  186.0,  188.0,  190.0,  192.0,  194.0,  196.0,  198.0,  200.0,
114 202.0,  204.0,  206.0,  208.0,  210.0,  212.0,  214.0,  216.0,  218.0,  220.0,
115 222.0,  224.0,  226.0,  228.0,  230.0,  232.0,  234.0,  236.0,  238.0,  240.0,
116 242.0,  243.0,  244.0,  245.0,  246.0,  247.0,  248.0,  248.5,  249.0,  249.5   } ;
117
118
119
120 AliTPCCorrection::AliTPCCorrection() 
121   : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
122 {
123   //
124   // default constructor
125   //
126 }
127
128 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
129 : TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
130 {
131   //
132   // default constructor, that set the name and title
133   //
134 }
135
136 AliTPCCorrection::~AliTPCCorrection() {
137   // 
138   // virtual destructor
139   //
140 }
141
142 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
143   //
144   // Corrects the initial coordinates x (cartesian coordinates)
145   // according to the given effect (inherited classes)
146   // roc represents the TPC read out chamber (offline numbering convention)
147   //
148   Float_t dx[3];
149   GetCorrection(x,roc,dx);
150   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
151 }
152
153 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
154   //
155   // Corrects the initial coordinates x (cartesian coordinates) and stores the new 
156   // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
157   // roc represents the TPC read out chamber (offline numbering convention)
158   //
159   Float_t dx[3];
160   GetCorrection(x,roc,dx);
161   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
162 }
163
164 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
165   //
166   // Distorts 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   GetDistortion(x,roc,dx);
172   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
173 }
174
175 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
176   //
177   // Distorts 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   GetDistortion(x,roc,dx);
183   for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
184 }
185
186 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
187   //
188   // This function delivers the correction values dx in respect to the inital coordinates x
189   // roc represents the TPC read out chamber (offline numbering convention)
190   // Note: The dx is overwritten by the inherited effectice class ...
191   //
192   for (Int_t j=0;j<3;++j) { dx[j]=0.; }
193 }
194
195 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
196   //
197   // This function delivers the distortion values dx in respect to the inital coordinates x
198   // roc represents the TPC read out chamber (offline numbering convention)
199   //
200   GetCorrection(x,roc,dx);
201   for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
202 }
203
204 void AliTPCCorrection::Init() {
205   //
206   // Initialization funtion (not used at the moment)
207   //
208 }
209
210 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
211   //
212   // Update function 
213   //
214 }
215
216 void AliTPCCorrection::Print(Option_t* /*option*/) const {
217   //
218   // Print function to check which correction classes are used 
219   // option=="d" prints details regarding the setted magnitude 
220   // option=="a" prints the C0 and C1 coefficents for calibration purposes
221   //
222   printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
223 }
224
225 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
226   //
227   // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
228   // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
229   // calibration run
230   //
231   fT1=t1;
232   fT2=t2;
233   //SetOmegaTauT1T2(omegaTau, t1, t2);
234 }
235
236 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
237   //
238   // Simple plot functionality.
239   // Returns a 2d hisogram which represents the corrections in radial direction (dr)
240   // in respect to position z within the XY plane.
241   // The histogramm has nx times ny entries. 
242   //
243   
244   TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
245                      nx,-250.,250.,ny,-250.,250.);
246   Float_t x[3],dx[3];
247   x[2]=z;
248   Int_t roc=z>0.?0:18; // FIXME
249   for (Int_t iy=1;iy<=ny;++iy) {
250     x[1]=h->GetYaxis()->GetBinCenter(iy);
251     for (Int_t ix=1;ix<=nx;++ix) {
252       x[0]=h->GetXaxis()->GetBinCenter(ix);
253       GetCorrection(x,roc,dx);
254       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
255       if (90.<=r0 && r0<=250.) {
256         Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
257         h->SetBinContent(ix,iy,r1-r0);
258       }
259       else
260         h->SetBinContent(ix,iy,0.);
261     }
262   }
263   return h;
264 }
265
266 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
267   //
268   // Simple plot functionality.
269   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi) 
270   // in respect to position z within the XY plane.
271   // The histogramm has nx times ny entries. 
272   //
273
274   TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
275                      nx,-250.,250.,ny,-250.,250.);
276   Float_t x[3],dx[3];
277   x[2]=z;
278   Int_t roc=z>0.?0:18; // FIXME
279   for (Int_t iy=1;iy<=ny;++iy) {
280     x[1]=h->GetYaxis()->GetBinCenter(iy);
281     for (Int_t ix=1;ix<=nx;++ix) {
282       x[0]=h->GetXaxis()->GetBinCenter(ix);
283       GetCorrection(x,roc,dx);
284       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
285       if (90.<=r0 && r0<=250.) {
286         Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
287         Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
288
289         Float_t dphi=phi1-phi0;
290         if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
291         if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
292       
293         h->SetBinContent(ix,iy,r0*dphi);
294       }
295       else
296         h->SetBinContent(ix,iy,0.);
297     }
298   }
299   return h;
300 }
301
302 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
303   //
304   // Simple plot functionality.
305   // Returns a 2d hisogram which represents the corrections in r direction (dr) 
306   // in respect to angle phi within the ZR plane.
307   // The histogramm has nx times ny entries. 
308   //
309   TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
310                      nz,-250.,250.,nr,85.,250.);
311   Float_t x[3],dx[3];
312   for (Int_t ir=1;ir<=nr;++ir) {
313     Float_t radius=h->GetYaxis()->GetBinCenter(ir);
314     x[0]=radius*TMath::Cos(phi);
315     x[1]=radius*TMath::Sin(phi);
316     for (Int_t iz=1;iz<=nz;++iz) {
317       x[2]=h->GetXaxis()->GetBinCenter(iz);
318       Int_t roc=x[2]>0.?0:18; // FIXME
319       GetCorrection(x,roc,dx);
320       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
321       Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
322       h->SetBinContent(iz,ir,r1-r0);
323     }
324   }
325   printf("SDF\n");
326   return h;
327
328 }
329
330 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
331   //
332   // Simple plot functionality.
333   // Returns a 2d hisogram which represents the corrections in rphi direction (drphi) 
334   // in respect to angle phi within the ZR plane.
335   // The histogramm has nx times ny entries. 
336   //
337   TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
338                      nz,-250.,250.,nr,85.,250.);
339   Float_t x[3],dx[3];
340   for (Int_t iz=1;iz<=nz;++iz) {
341     x[2]=h->GetXaxis()->GetBinCenter(iz);
342     Int_t roc=x[2]>0.?0:18; // FIXME
343     for (Int_t ir=1;ir<=nr;++ir) {
344       Float_t radius=h->GetYaxis()->GetBinCenter(ir);
345       x[0]=radius*TMath::Cos(phi);
346       x[1]=radius*TMath::Sin(phi);
347       GetCorrection(x,roc,dx);
348       Float_t r0=TMath::Sqrt((x[0]      )*(x[0]      )+(x[1]      )*(x[1]      ));
349       Float_t phi0=TMath::ATan2(x[1]      ,x[0]      );
350       Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
351       
352       Float_t dphi=phi1-phi0;
353       if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
354       if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
355       
356       h->SetBinContent(iz,ir,r0*dphi);
357     }
358   }
359   return h;
360 }
361
362 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
363                                    const char *xlabel,const char *ylabel,const char *zlabel,
364                                   Int_t nbinsx,Double_t xlow,Double_t xup,
365                                   Int_t nbinsy,Double_t ylow,Double_t yup) {
366   //
367   // Helper function to create a 2d histogramm of given size
368   //
369   
370   TString hname=name;
371   Int_t i=0;
372   if (gDirectory) {
373     while (gDirectory->FindObject(hname.Data())) {
374       hname =name;
375       hname+="_";
376       hname+=i;
377       ++i;
378     }
379   }
380   TH2F *h=new TH2F(hname.Data(),title,
381                    nbinsx,xlow,xup,
382                    nbinsy,ylow,yup);
383   h->GetXaxis()->SetTitle(xlabel);
384   h->GetYaxis()->SetTitle(ylabel);
385   h->GetZaxis()->SetTitle(zlabel);
386   h->SetStats(0);
387   return h;
388 }
389
390
391 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
392
393 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z, 
394                                                   const Double_t er[kNZ][kNR], Double_t &erValue ) {
395   //
396   // Interpolate table - 2D interpolation
397   //
398   Double_t saveEr[10] ;
399
400   Search( kNZ,   fgkZList,  z,   fJLow   ) ;
401   Search( kNR,   fgkRList,  r,   fKLow   ) ;
402   if ( fJLow < 0 ) fJLow = 0 ;   // check if out of range
403   if ( fKLow < 0 ) fKLow = 0 ;
404   if ( fJLow + order  >=    kNZ - 1 ) fJLow =   kNZ - 1 - order ;
405   if ( fKLow + order  >=    kNR - 1 ) fKLow =   kNR - 1 - order ;
406
407   for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
408       saveEr[j-fJLow]     = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r )   ;
409   }
410   erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z )   ;
411
412 }
413
414
415 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[], 
416                                        const Int_t order, const Double_t x ) {
417   //
418   // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
419   //
420
421   Double_t y ;
422   if ( order == 2 ) {                // Quadratic Interpolation = 2 
423     y  = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ; 
424     y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ; 
425     y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ; 
426   } else {                           // Linear Interpolation = 1
427     y  = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
428   }
429
430   return (y);
431
432 }
433
434
435 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
436   //
437   // Search an ordered table by starting at the most recently used point
438   //
439
440   Long_t middle, high ;
441   Int_t  ascend = 0, increment = 1 ;
442
443   if ( xArray[n-1] >= xArray[0] ) ascend = 1 ;  // Ascending ordered table if true
444   
445   if ( low < 0 || low > n-1 ) { 
446     low = -1 ; high = n ; 
447   } else {                                            // Ordered Search phase
448     if ( (Int_t)( x >= xArray[low] ) == ascend )  {
449       if ( low == n-1 ) return ;          
450       high = low + 1 ;
451       while ( (Int_t)( x >= xArray[high] ) == ascend ) {
452         low = high ;
453         increment *= 2 ;
454         high = low + increment ;
455         if ( high > n-1 )  {  high = n ; break ;  }
456       }
457     } else {
458       if ( low == 0 )  {  low = -1 ;  return ;  }
459       high = low - 1 ;
460       while ( (Int_t)( x < xArray[low] ) == ascend ) {
461         high = low ;
462         increment *= 2 ;
463         if ( increment >= high )  {  low = -1 ;  break ;  }
464         else  low = high - increment ;
465       }
466     }
467   }
468   
469   while ( (high-low) != 1 ) {                     // Binary Search Phase
470     middle = ( high + low ) / 2 ;
471     if ( (Int_t)( x >= xArray[middle] ) == ascend )
472       low = middle ;
473     else
474       high = middle ;
475   }
476   
477   if ( x == xArray[n-1] ) low = n-2 ;
478   if ( x == xArray[0]   ) low = 0 ;
479   
480 }
481
482
483 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
484   //
485   // Fit the track parameters - without and with distortion
486   // 1. Space points in the TPC are simulated along the trajectory  
487   // 2. Space points distorted
488   // 3. Fits  the non distorted and distroted track to the reference plane at refX
489   // 4. For visualization and debugging  purposes the space points and tracks can be stored  in the tree - using the TTreeSRedirector functionality 
490   //
491   // trackIn   - input track parameters
492   // refX     - reference X to fit the track
493   // dir      - direction - out=1 or in=-1
494   // pcstream -  debug streamer to check the results
495   //
496   // see AliExternalTrackParam.h documentation:
497   // track1.fP[0] - local y (rphi)
498   // track1.fP[1] - z 
499   // track1.fP[2] - sinus of local inclination angle
500   // track1.fP[3] - tangent of deep angle
501   // track1.fP[4] - 1/pt
502   AliTPCROC * roc = AliTPCROC::Instance();
503   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
504   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
505   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
506     
507   const Double_t kMaxSnp = 0.85;  
508   const Double_t kSigmaY=0.1;
509   const Double_t kSigmaZ=0.1;
510   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
511
512   AliExternalTrackParam  track(trackIn); // 
513   // generate points
514   AliTrackPointArray pointArray0(npoints0);
515   AliTrackPointArray pointArray1(npoints0);
516   Double_t xyz[3];
517   AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
518   //
519   // simulate the track
520   Int_t npoints=0;
521   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
522   for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
523     AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp);
524     track.GetXYZ(xyz);
525     xyz[0]+=gRandom->Gaus(0,0.001);
526     xyz[1]+=gRandom->Gaus(0,0.001);
527     AliTrackPoint pIn0;                               // space point          
528     AliTrackPoint pIn1;
529     Int_t sector= (xyz[2]>0)? 0:18;
530     pointArray0.GetPoint(pIn0,npoints);
531     pointArray1.GetPoint(pIn1,npoints);
532     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
533     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
534     DistortPoint(distPoint, sector);
535     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
536     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
537     //
538     track.Rotate(alpha);
539     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
540     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
541     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
542     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
543     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
544     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
545     pointArray0.AddPoint(npoints, &pIn0);      
546     pointArray1.AddPoint(npoints, &pIn1);
547     npoints++;
548     if (npoints>=npoints0) break;
549   }
550   //
551   // refit track
552   //
553   AliExternalTrackParam *track0=0;
554   AliExternalTrackParam *track1=0;
555   AliTrackPoint   point1,point2,point3;
556   if (dir==1) {  //make seed inner
557     pointArray0.GetPoint(point1,1);
558     pointArray0.GetPoint(point2,30);
559     pointArray0.GetPoint(point3,60);
560   }
561   if (dir==-1){ //make seed outer
562     pointArray0.GetPoint(point1,npoints-60);
563     pointArray0.GetPoint(point2,npoints-30);
564     pointArray0.GetPoint(point3,npoints-1);
565   }  
566   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
567   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
568
569
570   for (Int_t jpoint=0; jpoint<npoints; jpoint++){
571     Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
572     //
573     AliTrackPoint pIn0;
574     AliTrackPoint pIn1;
575     pointArray0.GetPoint(pIn0,ipoint);
576     pointArray1.GetPoint(pIn1,ipoint);
577     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
578     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
579     //
580     AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
581     AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
582     track.GetXYZ(xyz);  // distorted track also propagated to the same reference radius
583     //
584     Double_t pointPos[2]={0,0};
585     Double_t pointCov[3]={0,0,0};
586     pointPos[0]=prot0.GetY();//local y
587     pointPos[1]=prot0.GetZ();//local z
588     pointCov[0]=prot0.GetCov()[3];//simay^2
589     pointCov[1]=prot0.GetCov()[4];//sigmayz
590     pointCov[2]=prot0.GetCov()[5];//sigmaz^2
591     track0->Update(pointPos,pointCov);
592     //
593     Double_t deltaX=prot1.GetX()-prot0.GetX();   // delta X 
594     Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp()));  // deltaY due  delta X
595     Double_t deltaZX=deltaX*track1->GetTgl();                           // deltaZ due  delta X
596
597     pointPos[0]=prot1.GetY()-deltaYX;//local y
598     pointPos[1]=prot1.GetZ()-deltaZX;//local z
599     pointCov[0]=prot1.GetCov()[3];//simay^2
600     pointCov[1]=prot1.GetCov()[4];//sigmayz
601     pointCov[2]=prot1.GetCov()[5];//sigmaz^2
602     track1->Update(pointPos,pointCov);
603   }
604
605   AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
606   track1->Rotate(track0->GetAlpha());
607   track1->PropagateTo(track0->GetX(),AliTrackerBase::GetBz());
608
609   if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
610     "point0.="<<&pointArray0<<   //  points
611     "point1.="<<&pointArray1<<   //  distorted points
612     "trackIn.="<<&track<<       //  original track
613     "track0.="<<track0<<         //  fitted track
614     "track1.="<<track1<<         //  fitted distorted track
615     "\n";
616   new(&trackIn) AliExternalTrackParam(*track0);
617   delete track0;
618   return track1;
619 }
620
621
622
623
624
625 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
626   //
627   // create the distortion tree on a mesh with granularity given by step
628   // return the tree with distortions at given position 
629   // Map is created on the mesh with given step size
630   //
631   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
632   Float_t xyz[3];
633   for (Double_t x= -250; x<250; x+=step){
634     for (Double_t y= -250; y<250; y+=step){
635       Double_t r    = TMath::Sqrt(x*x+y*y);
636       if (r<80) continue;
637       if (r>250) continue;      
638       for (Double_t z= -250; z<250; z+=step){
639         Int_t roc=(z>0)?0:18;
640         xyz[0]=x;
641         xyz[1]=y;
642         xyz[2]=z;
643         Double_t phi  = TMath::ATan2(y,x);
644         DistortPoint(xyz,roc);
645         Double_t r1    = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
646         Double_t phi1  = TMath::ATan2(xyz[1],xyz[0]);
647         if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
648         if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
649         Double_t dx = xyz[0]-x;
650         Double_t dy = xyz[1]-y;
651         Double_t dz = xyz[2]-z;
652         Double_t dr=r1-r;
653         Double_t drphi=(phi1-phi)*r;
654         (*pcstream)<<"distortion"<<
655           "x="<<x<<           // original position        
656           "y="<<y<<
657           "z="<<z<<
658           "r="<<r<<
659           "phi="<<phi<<   
660           "x1="<<xyz[0]<<      // distorted position
661           "y1="<<xyz[1]<<
662           "z1="<<xyz[2]<<
663           "r1="<<r1<<
664           "phi1="<<phi1<<
665           //
666           "dx="<<dx<<          // delta position
667           "dy="<<dy<<
668           "dz="<<dz<<
669           "dr="<<dr<<
670           "drphi="<<drphi<<
671           "\n";
672       }
673     }   
674   }
675   delete pcstream;
676   TFile f(Form("correction%s.root",GetName()));
677   TTree * tree = (TTree*)f.Get("distortion");
678   TTree * tree2= tree->CopyTree("1");
679   tree2->SetName(Form("dist%s",GetName()));
680   tree2->SetDirectory(0);
681   delete tree;
682   return tree2;
683 }
684
685
686
687
688 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
689   //
690   // Make a fit tree:
691   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
692   // calculates partial distortions
693   // Partial distortion is stored in the resulting tree
694   // Output is storred in the file distortion_<dettype>_<partype>.root
695   // Partial  distortion is stored with the name given by correction name
696   //
697   //
698   // Parameters of function:
699   // input     - input tree
700   // dtype     - distortion type 0 - ITSTPC,  1 -TPCTRD, 2 - TPCvertex 
701   // ppype     - parameter type
702   // corrArray - array with partial corrections
703   // step      - skipe entries  - if 1 all entries processed - it is slow
704   // debug     0 if debug on also space points dumped - it is slow
705   const Double_t kMaxSnp = 0.85;  
706   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
707   //  const Double_t kB2C=-0.299792458e-3;
708   const Int_t kMinEntries=50;
709   Double_t phi,theta, snp, mean,rms, entries;
710   tinput->SetBranchAddress("theta",&theta);
711   tinput->SetBranchAddress("phi", &phi);
712   tinput->SetBranchAddress("snp",&snp);
713   tinput->SetBranchAddress("mean",&mean);
714   tinput->SetBranchAddress("rms",&rms);
715   tinput->SetBranchAddress("entries",&entries);
716   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
717   //
718   Int_t nentries=tinput->GetEntries();
719   Int_t ncorr=corrArray->GetEntries();
720   Double_t corrections[100]; //
721   Double_t tPar[5];
722   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
723   Double_t refX=0;
724   Int_t dir=0;
725   if (dtype==0) {refX=85; dir=-1;}
726   if (dtype==1) {refX=275; dir=1;}
727   if (dtype==2) {refX=85; dir=-1;}
728   //
729   for (Int_t ientry=0; ientry<nentries; ientry+=step){
730     tinput->GetEntry(ientry);
731     tPar[0]=0;
732     tPar[1]=theta*refX;
733     tPar[2]=snp;
734     tPar[3]=theta;
735     tPar[4]=(gRandom->Rndm()-0.5)*0.02;  // should be calculated - non equal to 0
736     Double_t bz=AliTrackerBase::GetBz();
737     if (refX>10. && TMath::Abs(bz)>0.1 )  tPar[4]=snp/(refX*bz*kB2C*2);
738     tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
739     if (TMath::Abs(snp)>0.251) continue;
740     (*pcstream)<<"fit"<<
741       "bz="<<bz<<         // magnetic filed used
742       "dtype="<<dtype<<   // detector match type
743       "ptype="<<ptype<<   // parameter type
744       "theta="<<theta<<   // theta
745       "phi="<<phi<<       // phi 
746       "snp="<<snp<<       // snp
747       "mean="<<mean<<     // mean dist value
748       "rms="<<rms<<       // rms
749       "entries="<<entries;// number of entries in bin
750     //
751     for (Int_t icorr=0; icorr<ncorr; icorr++) {
752       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
753       corrections[icorr]=0;
754       if (entries>kMinEntries){
755         AliExternalTrackParam trackIn(refX,phi,tPar,cov);
756         AliExternalTrackParam *trackOut = 0;
757         if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
758         if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
759         if (dtype==0) {refX=85; dir=-1;}
760         if (dtype==1) {refX=275; dir=1;}
761         if (dtype==2) {refX=0; dir=-1;}
762         //
763         AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
764         AliTrackerBase::PropagateTrackToBxByBz(trackOut,refX,kMass,3,kTRUE,kMaxSnp);
765         //
766         corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
767         delete trackOut;      
768       }      
769       (*pcstream)<<"fit"<<
770         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
771     }
772     (*pcstream)<<"fit"<<"\n";
773   }
774   delete pcstream;
775 }
776
777
778
779
780 void   AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
781   //
782   // make a distortion map out ou fthe residual histogram
783   // Results are written to the debug streamer - pcstream
784   // Parameters:
785   //   his0       - input (4D) residual histogram
786   //   pcstream   - file to write the tree
787   //   run        - run number
788   // marian.ivanov@cern.ch
789   const Int_t kMinEntries=50;
790   Int_t nbins1=his0->GetAxis(1)->GetNbins();
791   Int_t first1=his0->GetAxis(1)->GetFirst();
792   Int_t last1 =his0->GetAxis(1)->GetLast();
793   //
794   Double_t bz=AliTrackerBase::GetBz();
795   Int_t idim[4]={0,1,2,3};
796   for (Int_t ibin1=first1; ibin1<last1; ibin1++){   //axis 1 - theta
797     //
798     his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
799     Double_t       x1= his0->GetAxis(1)->GetBinCenter(ibin1);
800     THnSparse * his1 = his0->Projection(4,idim);  // projected histogram according range1
801     Int_t nbins3     = his1->GetAxis(3)->GetNbins();
802     Int_t first3     = his1->GetAxis(3)->GetFirst();
803     Int_t last3      = his1->GetAxis(3)->GetLast();
804     //
805
806     for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){   // axis 3 - local angle
807       his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
808       Double_t      x3= his1->GetAxis(3)->GetBinCenter(ibin3);
809       if (ibin3<first3) {
810         his1->GetAxis(3)->SetRangeUser(-1,1);
811         x3=0;
812       }
813       THnSparse * his3= his1->Projection(4,idim);         //projected histogram according selection 3
814       Int_t nbins2    = his3->GetAxis(2)->GetNbins();
815       Int_t first2    = his3->GetAxis(2)->GetFirst();
816       Int_t last2     = his3->GetAxis(2)->GetLast();
817       //
818       for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
819         his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
820         Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
821         TH1 * hisDelta = his3->Projection(0);
822         //
823         Double_t entries = hisDelta->GetEntries();
824         Double_t mean=0, rms=0;
825         if (entries>kMinEntries){
826           mean    = hisDelta->GetMean(); 
827           rms = hisDelta->GetRMS(); 
828         }
829         (*pcstream)<<hname<<
830           "run="<<run<<
831           "bz="<<bz<<
832           "theta="<<x1<<
833           "phi="<<x2<<
834           "snp="<<x3<<
835           "entries="<<entries<<
836           "mean="<<mean<<
837           "rms="<<rms<<
838           "\n";
839         delete hisDelta;
840         printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
841       }
842       delete his3;
843     }
844     delete his1;
845   }
846 }
847
848
849
850
851
852 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
853   //
854   // Store object in the OCDB
855   // By default the object is stored in the current directory 
856   // default comment consit of user name and the date
857   //
858   TString ocdbStorage="";
859   ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
860   AliCDBMetaData *metaData= new AliCDBMetaData();
861   metaData->SetObjectClassName("AliTPCCorrection");
862   metaData->SetResponsible("Marian Ivanov");
863   metaData->SetBeamPeriod(1);
864   metaData->SetAliRootVersion("05-25-01"); //root version
865   TString userName=gSystem->GetFromPipe("echo $USER");
866   TString date=gSystem->GetFromPipe("date");
867
868   if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
869   if (comment) metaData->SetComment(comment);
870   AliCDBId* id1=NULL;
871   id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
872   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
873   gStorage->Put(this, (*id1), metaData);
874 }
875