]> 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
AliMixInputHandlerInfo : BeginEvent function for each input handler after fChain...
[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
064244d8 31#include "AliTPCcalibDB.h"
6c104224 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 //
064244d8 232 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
233 if (!clparam) {
234 AliFatal("TPC OCDB not initialized");
235 return 0;
236 }
15e48021 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;
6c104224 269}
270
15e48021 271
272
273
6c104224 274Double_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;
15e48021 295 if (dout<0) return 0;
6c104224 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
15e48021 318Double_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);
6acb8ddc 339 if (dedge-params[1]<0) return 0;
15e48021 340 if (dedge>10.*params[1]) return 1;
6acb8ddc 341 result = 1.-TMath::Exp(-params[0]*(dedge-params[1]));
15e48021 342 return result;
343}
344
345Double_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
349Double_t AliTPCPointCorrection::SGetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
350 //
351 //
352 return fgInstance->GetEdgeQ0(sector, padrow, y);
353}
354
355void 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));
774a5ee9 361 TMatrixD * mat1Covar = (TMatrixD*)(sideACov.At(isec));
15e48021 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
774a5ee9 393Double_t AliTPCPointCorrection::GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/, Int_t quadrant){
15e48021 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
430Double_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
774a5ee9 440Double_t AliTPCPointCorrection::GetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/){
15e48021 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
463Double_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
6c104224 475
476