Corrected UInt_t <-> Int_t conversion
[u/mrichter/AliRoot.git] / TPC / AliTPCPointCorrection.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   Unlinearities fitting:
19
20   Model for Outer field cage:
21   Unlinearities at the edge aproximated using two exponential decays.
22  
23   dz = dz0(r,z) +dr(r,z)*tan(theta) 
24   dy =          +dr(r,z)*tan(phi)
25
26    
27   
28     
29 */
30
31 #include "AliTPCcalibDB.h"
32 #include "TLinearFitter.h"
33 #include "Riostream.h"
34 #include "TList.h"
35 #include "TMath.h"
36 #include "TCanvas.h"
37 #include "TFile.h"
38 #include "TF1.h"
39 #include "TVectorD.h"
40 #include "AliLog.h"
41 #include "AliTPCROC.h"
42 #include "AliTPCClusterParam.h"
43 #include "AliTPCPointCorrection.h"
44
45 AliTPCPointCorrection* AliTPCPointCorrection::fgInstance = 0;
46
47 ClassImp(AliTPCPointCorrection)
48
49 AliTPCPointCorrection::AliTPCPointCorrection():
50   TNamed(),
51   fParamsOutR(38),
52   fParamsOutZ(38),
53   fParamOutRVersion(0),
54   fErrorsOutR(38),
55   fErrorsOutZ(38),
56   fParamOutZVersion(0),
57   //
58   fXIO(0),
59   fXmiddle(0),
60   fXquadrant(0),
61   fArraySectorIntParam(36), // array of sector alignment parameters
62   fArraySectorIntCovar(36), // array of sector alignment covariances 
63   //
64   // Kalman filter for global alignment
65   //
66   fSectorParam(0),     // Kalman parameter   
67   fSectorCovar(0)     // Kalman covariance  
68
69 {
70   //
71   // Default constructor
72   //
73   AliTPCROC * roc = AliTPCROC::Instance();
74   fXquadrant = roc->GetPadRowRadii(36,53);
75   fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
76   fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
77 }
78
79 AliTPCPointCorrection::AliTPCPointCorrection(const Text_t *name, const Text_t *title):
80   TNamed(name,title),
81   fParamsOutR(38),
82   fParamsOutZ(38),
83   fParamOutRVersion(0),
84   fErrorsOutR(38),
85   fErrorsOutZ(38),
86   fParamOutZVersion(0),
87   //
88   // 
89   //
90   fXIO(0),
91   fXmiddle(0),
92   fXquadrant(0),
93   fArraySectorIntParam(36), // array of sector alignment parameters
94   fArraySectorIntCovar(36), // array of sector alignment covariances 
95   //
96   // Kalman filter for global alignment
97   //
98   fSectorParam(0),     // Kalman parameter   for A side
99   fSectorCovar(0)     // Kalman covariance  for A side 
100   
101 {
102   //
103   // Non default constructor
104   //
105   AliTPCROC * roc = AliTPCROC::Instance();
106   fXquadrant = roc->GetPadRowRadii(36,53);
107   fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
108   fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
109
110 }
111
112 AliTPCPointCorrection::~AliTPCPointCorrection(){
113   //
114   //
115   //
116 }
117
118
119 AliTPCPointCorrection* AliTPCPointCorrection::Instance()
120 {
121   //
122   // Singleton implementation
123   // Returns an instance of this class, it is created if neccessary
124   //
125   if (fgInstance == 0){
126     fgInstance = new AliTPCPointCorrection();
127   }
128   return fgInstance;
129 }
130
131
132
133 Double_t AliTPCPointCorrection::GetDrOut(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
134   //
135   //  return radial correction
136   //
137   if (fParamOutRVersion==0) return CorrectionOutR0(isGlobal, type,cx,cy,cz,sector);
138   return 0;
139 }
140
141 Double_t      AliTPCPointCorrection::SGetDrOut(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
142   //
143   // return radial correction - static function
144   // 
145   return fgInstance->GetDrOut(isGlobal, type,cx,cy,cz,sector);
146 }
147
148
149
150
151 Double_t AliTPCPointCorrection::GetDzOut(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
152   //
153   //
154   //
155   if (fParamOutZVersion==0) return CorrectionOutZ0(isGlobal, type,cx,cy,cz,sector);
156   return 0;
157 }
158
159
160 Double_t      AliTPCPointCorrection::SGetDzOut(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
161   //
162   //
163   //
164   return fgInstance->GetDzOut(isGlobal, type,cx,cy,cz,sector);
165 }
166
167
168
169
170 Double_t AliTPCPointCorrection::CorrectionOutR0(Bool_t isGlobal, Bool_t type,  Double_t cx, Double_t cy, Double_t cz, Int_t sector){
171   //
172   // return dR correction - for correction version 0 
173   // Parameters:
174   // isGlobal   - is the x in global frame
175   // type       - kTRUE - use sectors - kFALSE - only Side param
176   // cx, cy,cz  - cluster position
177   // sector     - sector number
178   if (isGlobal){
179     // recalculate sector if in global frame
180     Double_t alpha    = TMath::ATan2(cy,cx);
181     if (alpha<0)  alpha+=TMath::Pi()*2;
182     sector = Int_t(18*alpha/(TMath::Pi()*2));
183   }
184
185   if (type==kFALSE) sector=36+(sector%36>=18);  // side level parameters
186   // dout
187   Double_t radius = TMath::Sqrt(cx*cx+cy*cy);  
188   AliTPCROC * roc = AliTPCROC::Instance();
189   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
190   Double_t dout = xout-radius;  
191   if (dout<0) return 0;
192   //drift
193   Double_t dr   = 0.5 - TMath::Abs(cz/250.);
194   //
195   //
196   TVectorD * vec = GetParamOutR(sector);
197   if (!vec) return 0;
198   Double_t eout10 = TMath::Exp(-dout/10.);
199   Double_t eout5 = TMath::Exp(-dout/5.);                    
200   Double_t eoutp  = (eout10+eout5)*0.5;
201   Double_t eoutm  = (eout10-eout5)*0.5;
202   //
203   Double_t result=0;
204   result+=(*vec)[1]*eoutp;
205   result+=(*vec)[2]*eoutm;
206   result+=(*vec)[3]*eoutp*dr;
207   result+=(*vec)[4]*eoutm*dr;
208   result+=(*vec)[5]*eoutp*dr*dr;
209   result+=(*vec)[6]*eoutm*dr*dr;
210   return result;
211 }
212
213 Double_t AliTPCPointCorrection::RPhiCOGCorrection(Int_t sector, Int_t padrow, Float_t pad, Float_t cy, Float_t y, Float_t z, Float_t ky,Float_t qMax, Float_t threshold){
214   //
215   // Calculates COG corection in RPHI direction
216   // cluster and track position  y is supposed to be corrected before for other effects
217   // (e.g ExB and alignemnt)
218   // Rphi distortion dependeds on the distance to the edge-pad, distance to the wire edge and
219   // relative distance to the center of the pad. Therefore the y position is trnsfromed to the 
220   // pad coordiante frame (correction offset (ExB alignemnt) substracted). 
221   //   
222   // Input parameters:
223   //
224   // sector - sector number - 0-71  - cl.GetDetector()
225   // padrow - padrow number -0-63 -IROC 0-95 OROC - cl->GetRow()
226   // pad    - mean pad number  - cl->GetPad()
227   // cy     - cluster y        - cl->GetY()
228   //  y     - track -or cluster y
229   //  qMax  - cluster max charge - cl-.GetMax()
230   //  threshold - clusterer threshold
231   //
232   AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
233   if (!clparam) {
234     AliFatal("TPC OCDB not initialized"); 
235     return 0;
236   }
237   Int_t padtype=0;
238   if (sector>=36) padtype = (padrow>64)?2:1;
239   Double_t padwidth=(padtype==0)? 0.4:0.6;
240
241   Double_t sigma = clparam->GetRMS0(0,padtype,250-TMath::Abs(z),ky)/padwidth;
242   //
243   Int_t   npads   =  AliTPCROC::Instance()->GetNPads(sector,padrow);
244   Float_t yshift  =  TMath::Abs(cy)-TMath::Abs((-npads*0.5+pad)*padwidth);   // the clusters can be corrected before
245                                             // shift to undo correction
246   
247   y -= yshift*((y>0)?1.:-1.);                             // y in the sector coordinate
248   Double_t y0     = npads*0.5*padwidth;
249   Double_t dy     = (TMath::Abs(y0)-TMath::Abs(y))/padwidth-0.5;  // distance to  the center of pad0   
250   //
251   Double_t padcenter = TMath::Nint(dy);
252   Double_t sumw=0;
253   Double_t sumwi=0;
254   for (Double_t ip=padcenter-2.5; ip<=padcenter+2.5;ip++){
255     Double_t weightGaus = qMax*TMath::Exp(-(dy-ip)*(dy-ip)/(2*sigma*sigma));
256     Double_t ypad       = (ip-npads*0.5)*padwidth;
257     Double_t weightGain = GetEdgeQ0(sector,padrow,ypad);
258     //
259     Double_t weight = TMath::Nint(weightGaus*weightGain);
260     if (ip<0 &&weight> threshold) weight = threshold;  // this is done in cl finder
261     if (weight<0) continue;
262     sumw+=weight;
263     sumwi+=weight*(ip);    
264   }
265   Double_t result =0;
266   if (sumw>0) result = sumwi/sumw;
267   result = (result-dy)*padwidth;
268   return result;
269 }
270
271
272
273
274 Double_t AliTPCPointCorrection::CorrectionOutZ0(Bool_t isGlobal, Bool_t type,  Double_t cx, Double_t cy, Double_t cz, Int_t sector){
275   //
276   // return dR correction - for correction version 0 
277   // Parameters:
278   // isGlobal   - is the x in global frame
279   // type       - kTRUE - use sectors - kFALSE - only Side param
280   // cx, cy,cz  - cluster position
281   // sector     - sector number
282   if (isGlobal){
283     // recalculate sector if in global frame
284     Double_t alpha    = TMath::ATan2(cy,cx);
285     if (alpha<0)  alpha+=TMath::Pi()*2;
286     sector = Int_t(18*alpha/(TMath::Pi()*2));
287   }
288
289   if (type==kFALSE) sector=36+(sector%36>=18);  // side level parameters
290   // dout
291   Double_t radius = TMath::Sqrt(cx*cx+cy*cy);  
292   AliTPCROC * roc = AliTPCROC::Instance();
293   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
294   Double_t dout = xout-radius;
295   if (dout<0) return 0;
296   //drift
297   Double_t dr   = 0.5 - TMath::Abs(cz/250.);
298   //
299   //
300   TVectorD * vec = GetParamOutR(sector);
301   if (!vec) return 0;
302   Double_t eout10 = TMath::Exp(-dout/10.);
303   Double_t eout5 = TMath::Exp(-dout/5.);                    
304   Double_t eoutp  = (eout10+eout5)*0.5;
305   Double_t eoutm  = (eout10-eout5)*0.5;
306   //
307   Double_t result=0;
308   result+=(*vec)[1]*eoutp;
309   result+=(*vec)[2]*eoutm;
310   result+=(*vec)[3]*eoutp*dr;
311   result+=(*vec)[4]*eoutm*dr;
312   result+=(*vec)[5]*eoutp*dr*dr;
313   result+=(*vec)[6]*eoutm*dr*dr;
314   return result;
315
316 }
317
318 Double_t  AliTPCPointCorrection::GetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
319   //
320   /* TF1 fexp("fexp","1-exp(-[0]*(x-[1]))",0,20)
321           | param [0] | param [1]
322      IROC | 4.71413   | 1.39558
323      OROC1| 2.11437   | 1.52643
324      OROC2| 2.15082   | 1.53537 
325   */
326
327   Double_t params[2]={100,0};
328   if (sector<36){
329     params[0]=4.71413; params[1]=1.39558;
330   }
331   if (sector>=36){
332     if (padrow<64) {  params[0]=2.114; params[1]=1.526;}
333     if (padrow>=64){  params[0]=2.15; params[1]=1.535;}
334   }
335   Double_t result= 1;
336   Double_t xrow  = AliTPCROC::Instance()->GetPadRowRadii(sector,padrow);
337   Double_t ymax  = TMath::Tan(TMath::Pi()/18.)*xrow;
338   Double_t dedge = ymax-TMath::Abs(y);
339   if (dedge-params[1]<0)             return 0;
340   if (dedge>10.*params[1]) return 1;
341   result = 1.-TMath::Exp(-params[0]*(dedge-params[1]));
342   return result;
343 }
344
345 Double_t AliTPCPointCorrection::SRPhiCOGCorrection(Int_t sector, Int_t padrow, Float_t pad, Float_t cy, Float_t y, Float_t z, Float_t ky,Float_t qMax, Float_t threshold){
346   return fgInstance->RPhiCOGCorrection(sector, padrow, pad, cy, y,z, ky, qMax, threshold);
347 }
348
349 Double_t AliTPCPointCorrection::SGetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
350   //
351   //
352   return fgInstance->GetEdgeQ0(sector, padrow, y);
353 }
354
355 void     AliTPCPointCorrection::AddCorrectionSector(TObjArray & sideAPar, TObjArray &sideCPar, TObjArray & sideACov, TObjArray &sideCCov, Bool_t reset){
356   //
357   //
358   //
359   for (Int_t isec=0;isec<36;isec++){
360     TMatrixD * mat1 = (TMatrixD*)(sideAPar.At(isec));
361     TMatrixD * mat1Covar = (TMatrixD*)(sideACov.At(isec));
362     if (!mat1) continue;
363     TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
364     TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
365     if (!mat0) {
366       fArraySectorIntParam.AddAt(mat1->Clone(),isec); 
367       fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec); 
368       continue;
369     }
370     (*mat0Covar)=(*mat1Covar);      
371     if (reset)   (*mat0)=(*mat1);
372     if (!reset) (*mat0)+=(*mat1);
373   }
374
375   for (Int_t isec=0;isec<36;isec++){
376     TMatrixD * mat1 = (TMatrixD*)(sideCPar.At(isec));
377     TMatrixD * mat1Covar = (TMatrixD*)(sideCCov.At(isec));
378     if (!mat1) continue;
379     TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
380     TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
381     if (!mat0) {
382       fArraySectorIntParam.AddAt(mat1->Clone(),isec); 
383       fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec); 
384       continue;
385     }
386     (*mat0Covar)=(*mat1Covar);      
387     if (reset)   (*mat0)=(*mat1);
388     if (!reset) (*mat0)+=(*mat1);
389   }
390 }
391
392
393 Double_t AliTPCPointCorrection::GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/, Int_t quadrant){
394   //
395   // Get position correction for given sector
396   //
397
398   TMatrixD * param = (TMatrixD*)fArraySectorIntParam.At(sector%36);
399   if (!param) return 0;
400   if (quadrant<0){   //recalc quadrant if not specified
401     if (lx<fXIO) quadrant=0;
402     if(lx>fXIO){
403       if (lx<fXquadrant) {
404         if (ly<0) quadrant=1;
405         if (ly>0) quadrant=2;
406       }
407       if (lx>=fXquadrant) {
408         if (ly<0) quadrant=3;
409         if (ly>0) quadrant=4;
410       }
411     }    
412   }
413   Double_t a10 = (*param)(quadrant*6+0,0);
414   Double_t a20 = (*param)(quadrant*6+1,0);
415   Double_t a21 = (*param)(quadrant*6+2,0);
416   Double_t dx  = (*param)(quadrant*6+3,0);
417   Double_t dy  = (*param)(quadrant*6+4,0);
418   Double_t dz  = (*param)(quadrant*6+5,0);
419   Double_t deltaX = lx-fXmiddle;
420   if (coord==0) return dx;
421   if (coord==1) return dy+deltaX*a10;
422   if (coord==2) return dz+deltaX*a20+ly*a21;
423   if (coord==3) return a10;
424   if (coord==4) return a20;
425   if (coord==5) return a21;
426   //
427   return 0;
428 }
429
430 Double_t AliTPCPointCorrection::SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz, Int_t quadrant){
431   //
432   //
433   //
434   if (!Instance()) return 0;
435   return Instance()->GetCorrectionSector(coord,sector,lx,ly,lz, quadrant);
436 }
437
438
439
440 Double_t AliTPCPointCorrection::GetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/){
441   //
442   // Get position correction for given sector
443   //
444   if (!fSectorParam) return 0;
445   
446   Double_t a10 = (*fSectorParam)(sector*6+0,0);
447   Double_t a20 = (*fSectorParam)(sector*6+1,0);
448   Double_t a21 = (*fSectorParam)(sector*6+2,0);
449   Double_t dx  = (*fSectorParam)(sector*6+3,0);
450   Double_t dy  = (*fSectorParam)(sector*6+4,0);
451   Double_t dz  = (*fSectorParam)(sector*6+5,0);
452   Double_t deltaX = lx-fXIO;
453   //
454   if (coord==0) return dx;
455   if (coord==1) return dy+deltaX*a10;
456   if (coord==2) return dz+deltaX*a20+ly*a21;
457   if (coord==3) return a10;
458   if (coord==4) return a20;
459   if (coord==5) return a21;
460   return 0;
461 }
462
463 Double_t AliTPCPointCorrection::SGetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){
464   //
465   //
466   //
467   if (!Instance()) return 0;
468   return Instance()->GetCorrection(coord,sector,lx,ly,lz);
469 }
470
471
472
473
474
475
476