]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCPointCorrection.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New task to check TPC-ITS track prolongation eff with cosmics
[u/mrichter/AliRoot.git] / TPC / AliTPCPointCorrection.cxx
CommitLineData
6c104224 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
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"
15e48021 42#include "AliTPCClusterParam.h"
6c104224 43#include "AliTPCPointCorrection.h"
44
45AliTPCPointCorrection* AliTPCPointCorrection::fgInstance = 0;
46
47ClassImp(AliTPCPointCorrection)
48
49AliTPCPointCorrection::AliTPCPointCorrection():
50 TNamed(),
51 fParamsOutR(38),
52 fParamsOutZ(38),
53 fParamOutRVersion(0),
54 fErrorsOutR(38),
55 fErrorsOutZ(38),
15e48021 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
6c104224 69{
70 //
71 // Default constructor
72 //
15e48021 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;
6c104224 77}
78
79AliTPCPointCorrection::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),
15e48021 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
6c104224 101{
102 //
103 // Non default constructor
104 //
15e48021 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
6c104224 110}
111
112AliTPCPointCorrection::~AliTPCPointCorrection(){
113 //
114 //
115 //
116}
117
118
119AliTPCPointCorrection* 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
133Double_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
141Double_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
151Double_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
160Double_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
170Double_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);
15e48021 190 Double_t dout = xout-radius;
191 if (dout<0) return 0;
6c104224 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;
15e48021 211}
6c104224 212
15e48021 213Double_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 = AliTPCClusterParam::Instance();
233 Int_t padtype=0;
234 if (sector>=36) padtype = (padrow>64)?2:1;
235 Double_t padwidth=(padtype==0)? 0.4:0.6;
236
237 Double_t sigma = clparam->GetRMS0(0,padtype,250-TMath::Abs(z),ky)/padwidth;
238 //
239 Int_t npads = AliTPCROC::Instance()->GetNPads(sector,padrow);
240 Float_t yshift = TMath::Abs(cy)-TMath::Abs((-npads*0.5+pad)*padwidth); // the clusters can be corrected before
241 // shift to undo correction
242
243 y -= yshift*((y>0)?1.:-1.); // y in the sector coordinate
244 Double_t y0 = npads*0.5*padwidth;
245 Double_t dy = (TMath::Abs(y0)-TMath::Abs(y))/padwidth-0.5; // distance to the center of pad0
246 //
247 Double_t padcenter = TMath::Nint(dy);
248 Double_t sumw=0;
249 Double_t sumwi=0;
250 for (Double_t ip=padcenter-2.5; ip<=padcenter+2.5;ip++){
251 Double_t weightGaus = qMax*TMath::Exp(-(dy-ip)*(dy-ip)/(2*sigma*sigma));
252 Double_t ypad = (ip-npads*0.5)*padwidth;
253 Double_t weightGain = GetEdgeQ0(sector,padrow,ypad);
254 //
255 Double_t weight = TMath::Nint(weightGaus*weightGain);
256 if (ip<0 &&weight> threshold) weight = threshold; // this is done in cl finder
257 if (weight<0) continue;
258 sumw+=weight;
259 sumwi+=weight*(ip);
260 }
261 Double_t result =0;
262 if (sumw>0) result = sumwi/sumw;
263 result = (result-dy)*padwidth;
264 return result;
6c104224 265}
266
15e48021 267
268
269
6c104224 270Double_t AliTPCPointCorrection::CorrectionOutZ0(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
271 //
272 // return dR correction - for correction version 0
273 // Parameters:
274 // isGlobal - is the x in global frame
275 // type - kTRUE - use sectors - kFALSE - only Side param
276 // cx, cy,cz - cluster position
277 // sector - sector number
278 if (isGlobal){
279 // recalculate sector if in global frame
280 Double_t alpha = TMath::ATan2(cy,cx);
281 if (alpha<0) alpha+=TMath::Pi()*2;
282 sector = Int_t(18*alpha/(TMath::Pi()*2));
283 }
284
285 if (type==kFALSE) sector=36+(sector%36>=18); // side level parameters
286 // dout
287 Double_t radius = TMath::Sqrt(cx*cx+cy*cy);
288 AliTPCROC * roc = AliTPCROC::Instance();
289 Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
290 Double_t dout = xout-radius;
15e48021 291 if (dout<0) return 0;
6c104224 292 //drift
293 Double_t dr = 0.5 - TMath::Abs(cz/250.);
294 //
295 //
296 TVectorD * vec = GetParamOutR(sector);
297 if (!vec) return 0;
298 Double_t eout10 = TMath::Exp(-dout/10.);
299 Double_t eout5 = TMath::Exp(-dout/5.);
300 Double_t eoutp = (eout10+eout5)*0.5;
301 Double_t eoutm = (eout10-eout5)*0.5;
302 //
303 Double_t result=0;
304 result+=(*vec)[1]*eoutp;
305 result+=(*vec)[2]*eoutm;
306 result+=(*vec)[3]*eoutp*dr;
307 result+=(*vec)[4]*eoutm*dr;
308 result+=(*vec)[5]*eoutp*dr*dr;
309 result+=(*vec)[6]*eoutm*dr*dr;
310 return result;
311
312}
313
15e48021 314Double_t AliTPCPointCorrection::GetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
315 //
316 /* TF1 fexp("fexp","1-exp(-[0]*(x-[1]))",0,20)
317 | param [0] | param [1]
318 IROC | 4.71413 | 1.39558
319 OROC1| 2.11437 | 1.52643
320 OROC2| 2.15082 | 1.53537
321 */
322
323 Double_t params[2]={100,0};
324 if (sector<36){
325 params[0]=4.71413; params[1]=1.39558;
326 }
327 if (sector>=36){
328 if (padrow<64) { params[0]=2.114; params[1]=1.526;}
329 if (padrow>=64){ params[0]=2.15; params[1]=1.535;}
330 }
331 Double_t result= 1;
332 Double_t xrow = AliTPCROC::Instance()->GetPadRowRadii(sector,padrow);
333 Double_t ymax = TMath::Tan(TMath::Pi()/18.)*xrow;
334 Double_t dedge = ymax-TMath::Abs(y);
335 if (dedge-params[2]<0) return 0;
336 if (dedge>10.*params[1]) return 1;
337 result = 1.-TMath::Exp(-params[0]*(dedge-params[2]));
338 return result;
339}
340
341Double_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){
342 return fgInstance->RPhiCOGCorrection(sector, padrow, pad, cy, y,z, ky, qMax, threshold);
343}
344
345Double_t AliTPCPointCorrection::SGetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
346 //
347 //
348 return fgInstance->GetEdgeQ0(sector, padrow, y);
349}
350
351void AliTPCPointCorrection::AddCorrectionSector(TObjArray & sideAPar, TObjArray &sideCPar, TObjArray & sideACov, TObjArray &sideCCov, Bool_t reset){
352 //
353 //
354 //
355 for (Int_t isec=0;isec<36;isec++){
356 TMatrixD * mat1 = (TMatrixD*)(sideAPar.At(isec));
774a5ee9 357 TMatrixD * mat1Covar = (TMatrixD*)(sideACov.At(isec));
15e48021 358 if (!mat1) continue;
359 TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
360 TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
361 if (!mat0) {
362 fArraySectorIntParam.AddAt(mat1->Clone(),isec);
363 fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec);
364 continue;
365 }
366 (*mat0Covar)=(*mat1Covar);
367 if (reset) (*mat0)=(*mat1);
368 if (!reset) (*mat0)+=(*mat1);
369 }
370
371 for (Int_t isec=0;isec<36;isec++){
372 TMatrixD * mat1 = (TMatrixD*)(sideCPar.At(isec));
373 TMatrixD * mat1Covar = (TMatrixD*)(sideCCov.At(isec));
374 if (!mat1) continue;
375 TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
376 TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
377 if (!mat0) {
378 fArraySectorIntParam.AddAt(mat1->Clone(),isec);
379 fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec);
380 continue;
381 }
382 (*mat0Covar)=(*mat1Covar);
383 if (reset) (*mat0)=(*mat1);
384 if (!reset) (*mat0)+=(*mat1);
385 }
386}
387
388
774a5ee9 389Double_t AliTPCPointCorrection::GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/, Int_t quadrant){
15e48021 390 //
391 // Get position correction for given sector
392 //
393
394 TMatrixD * param = (TMatrixD*)fArraySectorIntParam.At(sector%36);
395 if (!param) return 0;
396 if (quadrant<0){ //recalc quadrant if not specified
397 if (lx<fXIO) quadrant=0;
398 if(lx>fXIO){
399 if (lx<fXquadrant) {
400 if (ly<0) quadrant=1;
401 if (ly>0) quadrant=2;
402 }
403 if (lx>=fXquadrant) {
404 if (ly<0) quadrant=3;
405 if (ly>0) quadrant=4;
406 }
407 }
408 }
409 Double_t a10 = (*param)(quadrant*6+0,0);
410 Double_t a20 = (*param)(quadrant*6+1,0);
411 Double_t a21 = (*param)(quadrant*6+2,0);
412 Double_t dx = (*param)(quadrant*6+3,0);
413 Double_t dy = (*param)(quadrant*6+4,0);
414 Double_t dz = (*param)(quadrant*6+5,0);
415 Double_t deltaX = lx-fXmiddle;
416 if (coord==0) return dx;
417 if (coord==1) return dy+deltaX*a10;
418 if (coord==2) return dz+deltaX*a20+ly*a21;
419 if (coord==3) return a10;
420 if (coord==4) return a20;
421 if (coord==5) return a21;
422 //
423 return 0;
424}
425
426Double_t AliTPCPointCorrection::SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz, Int_t quadrant){
427 //
428 //
429 //
430 if (!Instance()) return 0;
431 return Instance()->GetCorrectionSector(coord,sector,lx,ly,lz, quadrant);
432}
433
434
435
774a5ee9 436Double_t AliTPCPointCorrection::GetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/){
15e48021 437 //
438 // Get position correction for given sector
439 //
440 if (!fSectorParam) return 0;
441
442 Double_t a10 = (*fSectorParam)(sector*6+0,0);
443 Double_t a20 = (*fSectorParam)(sector*6+1,0);
444 Double_t a21 = (*fSectorParam)(sector*6+2,0);
445 Double_t dx = (*fSectorParam)(sector*6+3,0);
446 Double_t dy = (*fSectorParam)(sector*6+4,0);
447 Double_t dz = (*fSectorParam)(sector*6+5,0);
448 Double_t deltaX = lx-fXIO;
449 //
450 if (coord==0) return dx;
451 if (coord==1) return dy+deltaX*a10;
452 if (coord==2) return dz+deltaX*a20+ly*a21;
453 if (coord==3) return a10;
454 if (coord==4) return a20;
455 if (coord==5) return a21;
456 return 0;
457}
458
459Double_t AliTPCPointCorrection::SGetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){
460 //
461 //
462 //
463 if (!Instance()) return 0;
464 return Instance()->GetCorrection(coord,sector,lx,ly,lz);
465}
466
467
468
469
470
6c104224 471
472