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