1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Unlinearities fitting:
20 Model for Outer field cage:
21 Unlinearities at the edge aproximated using two exponential decays.
23 dz = dz0(r,z) +dr(r,z)*tan(theta)
24 dy = +dr(r,z)*tan(phi)
31 #include "AliTPCcalibDB.h"
32 #include "TLinearFitter.h"
33 #include "Riostream.h"
41 #include "AliTPCROC.h"
42 #include "AliTPCClusterParam.h"
43 #include "AliTPCPointCorrection.h"
45 AliTPCPointCorrection* AliTPCPointCorrection::fgInstance = 0;
47 ClassImp(AliTPCPointCorrection)
49 AliTPCPointCorrection::AliTPCPointCorrection():
61 fArraySectorIntParam(36), // array of sector alignment parameters
62 fArraySectorIntCovar(36), // array of sector alignment covariances
64 // Kalman filter for global alignment
66 fSectorParam(0), // Kalman parameter
67 fSectorCovar(0) // Kalman covariance
71 // Default constructor
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;
79 AliTPCPointCorrection::AliTPCPointCorrection(const Text_t *name, const Text_t *title):
93 fArraySectorIntParam(36), // array of sector alignment parameters
94 fArraySectorIntCovar(36), // array of sector alignment covariances
96 // Kalman filter for global alignment
98 fSectorParam(0), // Kalman parameter for A side
99 fSectorCovar(0) // Kalman covariance for A side
103 // Non default constructor
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;
112 AliTPCPointCorrection::~AliTPCPointCorrection(){
119 AliTPCPointCorrection* AliTPCPointCorrection::Instance()
122 // Singleton implementation
123 // Returns an instance of this class, it is created if neccessary
125 if (fgInstance == 0){
126 fgInstance = new AliTPCPointCorrection();
133 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
137 if (fParamOutRVersion==0) return CorrectionOutR0(isGlobal, type,cx,cy,cz,sector);
141 Double_t AliTPCPointCorrection::SGetDrOut(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
143 // return radial correction - static function
145 return fgInstance->GetDrOut(isGlobal, type,cx,cy,cz,sector);
151 Double_t AliTPCPointCorrection::GetDzOut(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
155 if (fParamOutZVersion==0) return CorrectionOutZ0(isGlobal, type,cx,cy,cz,sector);
160 Double_t AliTPCPointCorrection::SGetDzOut(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
164 return fgInstance->GetDzOut(isGlobal, type,cx,cy,cz,sector);
170 Double_t AliTPCPointCorrection::CorrectionOutR0(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
172 // return dR correction - for correction version 0
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
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));
185 if (type==kFALSE) sector=36+(sector%36>=18); // side level parameters
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;
193 Double_t dr = 0.5 - TMath::Abs(cz/250.);
196 TVectorD * vec = GetParamOutR(sector);
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;
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;
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){
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).
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
232 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
234 AliFatal("TPC OCDB not initialized");
238 if (sector>=36) padtype = (padrow>64)?2:1;
239 Double_t padwidth=(padtype==0)? 0.4:0.6;
241 Double_t sigma = clparam->GetRMS0(0,padtype,250-TMath::Abs(z),ky)/padwidth;
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
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
251 Double_t padcenter = TMath::Nint(dy);
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);
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;
266 if (sumw>0) result = sumwi/sumw;
267 result = (result-dy)*padwidth;
274 Double_t AliTPCPointCorrection::CorrectionOutZ0(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
276 // return dR correction - for correction version 0
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
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));
289 if (type==kFALSE) sector=36+(sector%36>=18); // side level parameters
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;
297 Double_t dr = 0.5 - TMath::Abs(cz/250.);
300 TVectorD * vec = GetParamOutR(sector);
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;
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;
318 Double_t AliTPCPointCorrection::GetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
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
327 Double_t params[2]={100,0};
329 params[0]=4.71413; params[1]=1.39558;
332 if (padrow<64) { params[0]=2.114; params[1]=1.526;}
333 if (padrow>=64){ params[0]=2.15; params[1]=1.535;}
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[2]<0) return 0;
340 if (dedge>10.*params[1]) return 1;
341 result = 1.-TMath::Exp(-params[0]*(dedge-params[2]));
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);
349 Double_t AliTPCPointCorrection::SGetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
352 return fgInstance->GetEdgeQ0(sector, padrow, y);
355 void AliTPCPointCorrection::AddCorrectionSector(TObjArray & sideAPar, TObjArray &sideCPar, TObjArray & sideACov, TObjArray &sideCCov, Bool_t reset){
359 for (Int_t isec=0;isec<36;isec++){
360 TMatrixD * mat1 = (TMatrixD*)(sideAPar.At(isec));
361 TMatrixD * mat1Covar = (TMatrixD*)(sideACov.At(isec));
363 TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
364 TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
366 fArraySectorIntParam.AddAt(mat1->Clone(),isec);
367 fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec);
370 (*mat0Covar)=(*mat1Covar);
371 if (reset) (*mat0)=(*mat1);
372 if (!reset) (*mat0)+=(*mat1);
375 for (Int_t isec=0;isec<36;isec++){
376 TMatrixD * mat1 = (TMatrixD*)(sideCPar.At(isec));
377 TMatrixD * mat1Covar = (TMatrixD*)(sideCCov.At(isec));
379 TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
380 TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
382 fArraySectorIntParam.AddAt(mat1->Clone(),isec);
383 fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec);
386 (*mat0Covar)=(*mat1Covar);
387 if (reset) (*mat0)=(*mat1);
388 if (!reset) (*mat0)+=(*mat1);
393 Double_t AliTPCPointCorrection::GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/, Int_t quadrant){
395 // Get position correction for given sector
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;
404 if (ly<0) quadrant=1;
405 if (ly>0) quadrant=2;
407 if (lx>=fXquadrant) {
408 if (ly<0) quadrant=3;
409 if (ly>0) quadrant=4;
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;
430 Double_t AliTPCPointCorrection::SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz, Int_t quadrant){
434 if (!Instance()) return 0;
435 return Instance()->GetCorrectionSector(coord,sector,lx,ly,lz, quadrant);
440 Double_t AliTPCPointCorrection::GetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/){
442 // Get position correction for given sector
444 if (!fSectorParam) return 0;
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;
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;
463 Double_t AliTPCPointCorrection::SGetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){
467 if (!Instance()) return 0;
468 return Instance()->GetCorrection(coord,sector,lx,ly,lz);