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