]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/Base/AliTPCCorrectionLookupTable.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
ATO-78 - edge handling (mirroring out of range pads) set as an option
[u/mrichter/AliRoot.git] / TPC / Base / AliTPCCorrectionLookupTable.cxx
CommitLineData
05da1b4e 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#include <TMath.h>
17#include <TMatrixF.h>
2d4e971f 18#include <TStopwatch.h>
8d79bbb7 19#include <TString.h>
20#include <TFile.h>
21#include <TObjArray.h>
22#include <TSystem.h>
23#include <THashList.h>
fd245273 24#include <TVector2.h>
25#include <TLinearFitter.h>
05da1b4e 26
27#include <AliLog.h>
28#include <AliTPCROC.h>
29
30#include "AliTPCCorrection.h"
31
32#include "AliTPCCorrectionLookupTable.h"
33
34ClassImp(AliTPCCorrectionLookupTable)
35
36//_________________________________________________________________________________________
37AliTPCCorrectionLookupTable::AliTPCCorrectionLookupTable()
38: AliTPCCorrection()
8d79bbb7 39, fNR(0)
40, fNPhi(0)
41, fNZ(0)
862220e2 42, fCorrScaleFactor(-1)
fd245273 43, fFillCorrection(kTRUE)
531a1f0b 44, fLimitsR()
45, fLimitsPhi()
46, fLimitsZ()
05da1b4e 47, fLookUpDxDist(0x0)
48, fLookUpDyDist(0x0)
49, fLookUpDzDist(0x0)
50, fLookUpDxCorr(0x0)
51, fLookUpDyCorr(0x0)
52, fLookUpDzCorr(0x0)
53{
54 //
55 //
56 //
57}
58//_________________________________________________________________________________________
59AliTPCCorrectionLookupTable::~AliTPCCorrectionLookupTable()
60{
61 //
62 // dtor
63 //
64
65 ResetTables();
66}
67
68//_________________________________________________________________________________________
69void AliTPCCorrectionLookupTable::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
70 //
71 // Get interpolated Distortion
72 //
73
2d4e971f 74 GetInterpolation(x,roc,dx,fLookUpDxDist,fLookUpDyDist,fLookUpDzDist);
05da1b4e 75}
76
77//_________________________________________________________________________________________
78void AliTPCCorrectionLookupTable::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
79 //
80 // Get interplolated correction
81 //
2d4e971f 82 GetInterpolation(x,roc,dx,fLookUpDxCorr,fLookUpDyCorr,fLookUpDzCorr);
862220e2 83
84 if (fCorrScaleFactor>0) {
85 dx[0]*=fCorrScaleFactor;
86 dx[1]*=fCorrScaleFactor;
87 dx[2]*=fCorrScaleFactor;
88 }
05da1b4e 89}
90
91//_________________________________________________________________________________________
92void AliTPCCorrectionLookupTable::GetInterpolation(const Float_t x[],const Short_t roc,Float_t dx[],
2d4e971f 93 TMatrixF **mDx, TMatrixF **mDy, TMatrixF **mDz)
05da1b4e 94{
95 //
96 // Calculates the correction/distotring from a lookup table
97 // type: 0 = correction
98 // 1 = distortion
99 //
100
101// Float_t typeSign=-1;
102// if (type==1) typeSign=1;
103
104 Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
105
106 Double_t r, phi, z ;
107 Int_t sign;
108
109 r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
110 phi = TMath::ATan2(x[1],x[0]) ;
111 if ( phi < 0 ) phi += TMath::TwoPi() ; // Table uses phi from 0 to 2*Pi
112 z = x[2] ; // Create temporary copy of x[2]
113
114 if ( (roc%36) < 18 ) {
115 sign = 1; // (TPC A side)
116 } else {
117 sign = -1; // (TPC C side)
118 }
119
120 if ( sign==1 && z < fgkZOffSet ) z = fgkZOffSet; // Protect against discontinuity at CE
121 if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet; // Protect against discontinuity at CE
122
123
124 if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
125 AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
126
127 // Get the Er and Ephi field integrals plus the integral over Z
128 dx[0] = Interpolate3DTable(order, r, z, phi,
129 fNR, fNZ, fNPhi,
531a1f0b 130 fLimitsR.GetMatrixArray(),
131 fLimitsZ.GetMatrixArray(),
132 fLimitsPhi.GetMatrixArray(),
05da1b4e 133 mDx );
134 dx[1] = Interpolate3DTable(order, r, z, phi,
135 fNR, fNZ, fNPhi,
531a1f0b 136 fLimitsR.GetMatrixArray(),
137 fLimitsZ.GetMatrixArray(),
138 fLimitsPhi.GetMatrixArray(),
05da1b4e 139 mDy);
140 dx[2] = Interpolate3DTable(order, r, z, phi,
141 fNR, fNZ, fNPhi,
531a1f0b 142 fLimitsR.GetMatrixArray(),
143 fLimitsZ.GetMatrixArray(),
144 fLimitsPhi.GetMatrixArray(),
05da1b4e 145 mDz );
146}
147
148//_________________________________________________________________________________________
2d4e971f 149void AliTPCCorrectionLookupTable::CreateLookupTable(AliTPCCorrection &tpcCorr, Float_t stepSize/*=5.*/)
05da1b4e 150{
151 //
8d79bbb7 152 // create lookup table for all phi,r,z bins
05da1b4e 153 //
154
8d79bbb7 155 if (fNR==0) {
156 AliError("Limits are not set yet. Please use one of the Set..Limits functions first");
157 return;
158 }
05da1b4e 159
2d4e971f 160 TStopwatch s;
161
162 ResetTables();
163 InitTables();
164
8d79bbb7 165 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
166 CreateLookupTablePhiBin(tpcCorr,iPhi,stepSize);
167 }
168
169 s.Stop();
170 AliInfo(Form("Required time for lookup table creation: %.2f (%.2f) sec. real (cpu)",s.RealTime(),s.CpuTime()));
171}
172
173//_________________________________________________________________________________________
174void AliTPCCorrectionLookupTable::CreateLookupTableSinglePhi(AliTPCCorrection &tpcCorr, Int_t iPhi, Float_t stepSize)
175{
176 //
177 // Lookup table for only one phi bin. Can be used for parallel processing
178 //
179
180 if (fNR==0) {
181 AliError("Limits are not set yet. Please use one of the Set..Limits functions first");
182 return;
183 }
184
185 TStopwatch s;
186
187 ResetTables();
188 InitTableArrays();
189 InitTablesPhiBin(iPhi);
190 CreateLookupTablePhiBin(tpcCorr,iPhi,stepSize);
191
192 s.Stop();
193 AliInfo(Form("Required time for lookup table creation: %.2f (%.2f) sec. real (cpu)",s.RealTime(),s.CpuTime()));
194}
195
196//_________________________________________________________________________________________
197void AliTPCCorrectionLookupTable::CreateLookupTablePhiBin(AliTPCCorrection &tpcCorr, Int_t iPhi, Float_t stepSize)
198{
199 //
200 //
201 //
202
203 if (iPhi<0||iPhi>=fNPhi) return;
204
2d4e971f 205 const Float_t delta=stepSize; // 5cm
05da1b4e 206 Float_t x[3]={0.,0.,0.};
207 Float_t dx[3]={0.,0.,0.};
8d79bbb7 208
209 Double_t phi=fLimitsPhi(iPhi);
210 //
211 TMatrixF &mDxDist = *fLookUpDxDist[iPhi];
212 TMatrixF &mDyDist = *fLookUpDyDist[iPhi];
213 TMatrixF &mDzDist = *fLookUpDzDist[iPhi];
214 //
215 TMatrixF &mDxCorr = *fLookUpDxCorr[iPhi];
216 TMatrixF &mDyCorr = *fLookUpDyCorr[iPhi];
217 TMatrixF &mDzCorr = *fLookUpDzCorr[iPhi];
05da1b4e 218
8d79bbb7 219 for (Int_t ir=0; ir<fNR; ++ir){
220 Double_t r=fLimitsR(ir);
221 x[0]=r * TMath::Cos(phi);
222 x[1]=r * TMath::Sin(phi);
05da1b4e 223
8d79bbb7 224 for (Int_t iz=0; iz<fNZ; ++iz){
225 Double_t z=fLimitsZ(iz);
226 x[2]=z;
227 //TODO: change hardcoded value for r>133.?
228 Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
229 if (r>133.) roc+=36;
230 if (z<0) roc+=18;
231
232 if (delta>0)
233 tpcCorr.GetDistortionIntegralDz(x,roc,dx,delta);
234 else
235 tpcCorr.GetDistortion(x,roc,dx);
236 mDxDist(ir,iz)=dx[0];
237 mDyDist(ir,iz)=dx[1];
238 mDzDist(ir,iz)=dx[2];
fd245273 239
240 if (fFillCorrection) {
241 if (delta>0)
242 tpcCorr.GetCorrectionIntegralDz(x,roc,dx,delta);
243 else
244 tpcCorr.GetCorrection(x,roc,dx);
245 mDxCorr(ir,iz)=dx[0];
246 mDyCorr(ir,iz)=dx[1];
247 mDzCorr(ir,iz)=dx[2];
248 }
05da1b4e 249 }
250 }
8d79bbb7 251
05da1b4e 252}
253
254//_________________________________________________________________________________________
255void AliTPCCorrectionLookupTable::InitTables()
8d79bbb7 256{
257 //
258 // Init all tables
259 //
260
261 InitTableArrays();
262 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
263 InitTablesPhiBin(iPhi);
264 }
265}
266
fef4baf6 267//_________________________________________________________________________________________
268void AliTPCCorrectionLookupTable::CreateLookupTableFromResidualDistortion(THn &resDist)
269{
270 //
271 // create lookup table from residual distortions stored in a 3d histogram
272 // assume dimensions are r, phi, z
273 //
274 if (fNR==0) {
275 AliError("Limits are not set yet. Please use one of the Set..Limits functions first");
276 return;
277 }
278
fef4baf6 279 ResetTables();
280 InitTables();
281
282 Double_t x[3]={0.,0.,0.};
283
284 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
285 const Double_t phi=fLimitsPhi(iPhi);
286 x[1]=phi;
287 //
288 TMatrixF &mDxDist = *fLookUpDxDist[iPhi];
289 TMatrixF &mDyDist = *fLookUpDyDist[iPhi];
290 TMatrixF &mDzDist = *fLookUpDzDist[iPhi];
291 //
292 TMatrixF &mDxCorr = *fLookUpDxCorr[iPhi];
293 TMatrixF &mDyCorr = *fLookUpDyCorr[iPhi];
294 TMatrixF &mDzCorr = *fLookUpDzCorr[iPhi];
295
296 for (Int_t ir=0; ir<fNR; ++ir){
297 const Double_t r=fLimitsR(ir);
298 x[0]=r;
299
300 for (Int_t iz=0; iz<fNZ; ++iz){
301 const Double_t z=fLimitsZ(iz);
302 x[2]=z;
303
304 const Double_t drphi = resDist.GetBinContent(resDist.GetBin(x));
305 Double_t dx[3]={0.,drphi,0.};
306
307 // transform rphi distortions (local y, so dy') to a global distortion
308 // assume no radial distortion (dx' = 0)
309 // assume no residual distortion in z for the moment
310 Double_t cs=TMath::Cos(phi), sn=TMath::Sin(phi), lx=dx[0];
311 dx[0]=lx*cs - dx[1]*sn; dx[1]=lx*sn + dx[1]*cs;
312
313 mDxDist(ir,iz)=dx[0];
314 mDyDist(ir,iz)=dx[1];
315 mDzDist(ir,iz)=dx[2];
316
317 mDxCorr(ir,iz)=-dx[0];
318 mDyCorr(ir,iz)=-dx[1];
319 mDzCorr(ir,iz)=-dx[2];
320 }
321 }
322 }
323}
637ba19a 324
325//_________________________________________________________________________________________
326void AliTPCCorrectionLookupTable::CreateResidual(AliTPCCorrection *distortion, AliTPCCorrection* correction)
327{
328 //
329 // create lookup table from residual distortions calculated from distorted - correction
330 //
331
332 ResetTables();
333 InitTables();
334
335 Float_t x[3]={0.,0.,0.};
fef4baf6 336
637ba19a 337 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
338 const Double_t phi=fLimitsPhi(iPhi);
339 //
340 TMatrixF &mDxDist = *fLookUpDxDist[iPhi];
341 TMatrixF &mDyDist = *fLookUpDyDist[iPhi];
342 TMatrixF &mDzDist = *fLookUpDzDist[iPhi];
343 //
344 TMatrixF &mDxCorr = *fLookUpDxCorr[iPhi];
345 TMatrixF &mDyCorr = *fLookUpDyCorr[iPhi];
346 TMatrixF &mDzCorr = *fLookUpDzCorr[iPhi];
347
348 for (Int_t ir=0; ir<fNR; ++ir){
349 const Double_t r=fLimitsR(ir);
350 x[0]=r * TMath::Cos(phi);
351 x[1]=r * TMath::Sin(phi);
352
353 for (Int_t iz=0; iz<fNZ; ++iz){
354 const Double_t z=fLimitsZ(iz);
355 x[2]=z;
356
357 //original point
d8fd14dd 358 Float_t xdc[3]={x[0], x[1], x[2]};
637ba19a 359
360 Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
361 if (r>133.) roc+=36;
362 if (z<0) roc+=18;
363
364 //get residual distortion
d8fd14dd 365 distortion->DistortPoint(xdc, roc);
366 correction->CorrectPoint(xdc, roc);
6e4c2f43 367 Float_t dx[3]={xdc[0]-x[0], xdc[1]-x[1], xdc[2]-x[2]};
637ba19a 368
369 mDxDist(ir,iz)=dx[0];
370 mDyDist(ir,iz)=dx[1];
371 mDzDist(ir,iz)=dx[2];
372
373 mDxCorr(ir,iz)=-dx[0];
374 mDyCorr(ir,iz)=-dx[1];
375 mDzCorr(ir,iz)=-dx[2];
376 }
377 }
378 }
379}
380
8d79bbb7 381//_________________________________________________________________________________________
382void AliTPCCorrectionLookupTable::InitTablesPhiBin(Int_t iPhi)
383{
384 //
385 //
386 //
387
388 // check if already initialised
389 if (iPhi<0||iPhi>=fNPhi) return;
390 if (fLookUpDxCorr[iPhi]) return;
391
392 fLookUpDxCorr[iPhi] = new TMatrixF(fNR,fNZ);
393 fLookUpDyCorr[iPhi] = new TMatrixF(fNR,fNZ);
394 fLookUpDzCorr[iPhi] = new TMatrixF(fNR,fNZ);
395
396 fLookUpDxDist[iPhi] = new TMatrixF(fNR,fNZ);
397 fLookUpDyDist[iPhi] = new TMatrixF(fNR,fNZ);
398 fLookUpDzDist[iPhi] = new TMatrixF(fNR,fNZ);
399
400}
401//_________________________________________________________________________________________
402void AliTPCCorrectionLookupTable::InitTableArrays()
05da1b4e 403{
404 //
405 //
406 //
8d79bbb7 407
408 // needs to be before the table creation to set the limits
409 SetupDefaultLimits();
410
05da1b4e 411 fLookUpDxCorr = new TMatrixF*[fNPhi];
412 fLookUpDyCorr = new TMatrixF*[fNPhi];
413 fLookUpDzCorr = new TMatrixF*[fNPhi];
414
415 fLookUpDxDist = new TMatrixF*[fNPhi];
416 fLookUpDyDist = new TMatrixF*[fNPhi];
417 fLookUpDzDist = new TMatrixF*[fNPhi];
418
419 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
8d79bbb7 420 fLookUpDxCorr[iPhi] = 0x0;
421 fLookUpDyCorr[iPhi] = 0x0;
422 fLookUpDzCorr[iPhi] = 0x0;
05da1b4e 423
8d79bbb7 424 fLookUpDxDist[iPhi] = 0x0;
425 fLookUpDyDist[iPhi] = 0x0;
426 fLookUpDzDist[iPhi] = 0x0;
05da1b4e 427 }
05da1b4e 428}
429
430//_________________________________________________________________________________________
431void AliTPCCorrectionLookupTable::ResetTables()
432{
433 //
434 // Reset the lookup tables
435 //
436
437 if (!fLookUpDxCorr) return;
438
439 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
440 delete fLookUpDxCorr[iPhi];
441 delete fLookUpDyCorr[iPhi];
442 delete fLookUpDzCorr[iPhi];
443
444 delete fLookUpDxDist[iPhi];
445 delete fLookUpDyDist[iPhi];
446 delete fLookUpDzDist[iPhi];
447 }
448
449 delete [] fLookUpDxCorr;
450 delete [] fLookUpDyCorr;
451 delete [] fLookUpDzCorr;
452
453 delete [] fLookUpDxDist;
454 delete [] fLookUpDyDist;
455 delete [] fLookUpDzDist;
456
457 fLookUpDxCorr = 0x0;
458 fLookUpDyCorr = 0x0;
459 fLookUpDzCorr = 0x0;
460
461 fLookUpDxDist = 0x0;
462 fLookUpDyDist = 0x0;
463 fLookUpDzDist = 0x0;
05da1b4e 464}
465
466//_________________________________________________________________________________________
467void AliTPCCorrectionLookupTable::SetupDefaultLimits()
468{
469 //
470 // Set default limits for tables
471 //
472
531a1f0b 473 fNR = kNR;
474 fNPhi = kNPhi;
475 fNZ = kNZ;
476 fLimitsR. ResizeTo(fNR);
477 fLimitsPhi.ResizeTo(fNPhi);
478 fLimitsZ. ResizeTo(fNZ);
479 fLimitsR. SetElements(fgkRList);
480 fLimitsPhi.SetElements(fgkPhiList);
481 fLimitsZ. SetElements(fgkZList);
05da1b4e 482}
8d79bbb7 483
484//_________________________________________________________________________________________
485void AliTPCCorrectionLookupTable::ResetLimits()
486{
487 fNR = 0;
488 fNPhi = 0;
489 fNZ = 0;
490
491 fLimitsR. ResizeTo(1);
492 fLimitsPhi.ResizeTo(1);
493 fLimitsZ. ResizeTo(1);
494}
495
496//_________________________________________________________________________________________
497void AliTPCCorrectionLookupTable::MergePhiTables(const char* files)
498{
499 //
500 // merge all lookup tables stored in 'files' with this one
501 // assume that each lookup table in each file has only one phi bin
502 //
503
504 ResetTables();
505 ResetLimits(); // use limits from the first file assuming they are all the same
506
507 TString sfiles=gSystem->GetFromPipe(Form("ls %s",files));
508 TObjArray *arrFiles=sfiles.Tokenize("\n");
509
510 for (Int_t i=0; i<arrFiles->GetEntriesFast();++i){
511 TFile f(arrFiles->At(i)->GetName());
512 if (!f.IsOpen() || f.IsZombie()) continue;
513 AliTPCCorrectionLookupTable *lt=dynamic_cast<AliTPCCorrectionLookupTable*>(f.Get(f.GetListOfKeys()->First()->GetName()));
514 if (!lt) {
515 AliError(Form("first object in file '%s' is not of type AliTPCCorrectionLookupTable!",f.GetName()));
516 continue;
517 }
518
519 if (!fNR) {
520 fNR = lt->fNR;
521 fNPhi = lt->fNPhi;
522 fNZ = lt->fNZ;
523 fLimitsR = lt->fLimitsR;
524 fLimitsZ = lt->fLimitsZ;
525 fLimitsPhi = lt->fLimitsPhi;
526 InitTableArrays();
527 } else {
528 if (fNR!=lt->fNR || fNPhi!=lt->fNPhi || fNZ!=lt->fNZ ){
529 AliError(Form("Current limits don't macht the ones in file '%s'",f.GetName()));
530 continue;
531 }
532 }
533
534 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi) {
535 if (!lt->fLookUpDxCorr[iPhi]) continue;
536
537 AliInfo(Form("Adding phi bin '%d' from file '%s'",iPhi,f.GetName()));
538
539 InitTablesPhiBin(iPhi);
540
541 *fLookUpDxDist[iPhi] = *lt->fLookUpDxDist[iPhi];
542 *fLookUpDyDist[iPhi] = *lt->fLookUpDyDist[iPhi];
543 *fLookUpDzDist[iPhi] = *lt->fLookUpDzDist[iPhi];
544 //
545 *fLookUpDxCorr[iPhi] = *lt->fLookUpDxCorr[iPhi];
546 *fLookUpDyCorr[iPhi] = *lt->fLookUpDyCorr[iPhi];
547 *fLookUpDzCorr[iPhi] = *lt->fLookUpDzCorr[iPhi];
548 break;
549 }
550 }
551
552 //check of all phi bins are initialised
553 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi) {
554 if (!fLookUpDxCorr[iPhi]) {
c0701530 555 AliFatal(Form("Phi bin '%d' not initialised from files!",iPhi));
8d79bbb7 556 }
557 }
558
559 delete arrFiles;
fd245273 560}
561
562//_________________________________________________________________________________________
563void AliTPCCorrectionLookupTable::BuildExactInverse()
564{
565 //
566 // this method build the exact inverse of the standard distortion map
567 // for the the distortion man first needs to be calculated
568 // then the correction map will be overwritten
569 //
570
571 Float_t x[3] = {0.,0.,0.};
572 Float_t x2[3] = {0.,0.,0.};
573 Float_t xref[3] = {0.,0.,0.};
574 Float_t xd[3] = {0.,0.,0.};
575 Float_t dx[3] = {0.,0.,0.};
576
577 // reset correction matrices
578 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
579 TMatrixF &mDxCorr = *fLookUpDxCorr[iPhi];
580 TMatrixF &mDyCorr = *fLookUpDyCorr[iPhi];
581 TMatrixF &mDzCorr = *fLookUpDzCorr[iPhi];
582
583 for (Int_t ir=0; ir<fNR; ++ir){
584 for (Int_t iz=0; iz<fNZ; ++iz){
585 mDxCorr(ir,iz) = -1000.;
586 mDyCorr(ir,iz) = -1000.;
587 mDzCorr(ir,iz) = -1000.;
588 }
589 }
590 }
591
592 // get interplolated corrections on standard grid
593 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
594 Double_t phi=fLimitsPhi(iPhi);
595 TMatrixF &mDxDist = *fLookUpDxDist[iPhi];
596 TMatrixF &mDyDist = *fLookUpDyDist[iPhi];
597 TMatrixF &mDzDist = *fLookUpDzDist[iPhi];
598
599 for (Int_t ir=0; ir<fNR; ++ir){
600 Double_t r=fLimitsR(ir);
601 x[0]=r * TMath::Cos(phi);
602 x[1]=r * TMath::Sin(phi);
603
604 for (Int_t iz=0; iz<fNZ; ++iz){
605 Double_t z=fLimitsZ(iz);
606 x[2]=z;
607
608 //TODO: change hardcoded value for r>133.?
609 Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
610 if (r>133.) roc+=36;
611 if (z<0) roc+=18;
612
613 dx[0] = mDxDist(ir,iz);
614 dx[1] = mDyDist(ir,iz);
615 dx[2] = mDzDist(ir,iz);
616
617 xd[0] = x[0]+dx[0];
618 xd[1] = x[1]+dx[1];
619 xd[2] = x[2]+dx[2];
620
621 const Double_t phid = TVector2::Phi_0_2pi(TMath::ATan2(xd[1],xd[0]));
622 const Double_t rd = TMath::Sqrt(xd[0]*xd[0] + xd[1]*xd[1]);
623 const Double_t zd = xd[2];
624
625 Int_t ilow = 0, jlow = 0, klow = 0 ;
626
627 Search( fLimitsR.GetNrows(), fLimitsR.GetMatrixArray(), rd, ilow ) ;
628 Search( fLimitsZ.GetNrows(), fLimitsZ.GetMatrixArray(), zd, jlow ) ;
629 Search( fLimitsPhi.GetNrows(), fLimitsPhi.GetMatrixArray(), phid, klow ) ;
630
631 if ( ilow < 0 ) ilow = 0 ; // check if out of range
632 if ( jlow < 0 ) jlow = 0 ;
633 if ( klow < 0 ) klow = 0 ;
634 if ( ilow >= fLimitsR.GetNrows()) ilow = fLimitsR.GetNrows() - 1;
635 if ( jlow >= fLimitsZ.GetNrows()) jlow = fLimitsZ.GetNrows() - 1;
636 if ( klow >= fLimitsPhi.GetNrows()) klow = fLimitsPhi.GetNrows() - 1;
637
638 const Double_t phiRef = fLimitsPhi[klow];
639 const Double_t rRef = fLimitsR[ilow];
640 const Double_t zRef = fLimitsZ[jlow];
641
642 TMatrixF &mDxCorr = *fLookUpDxCorr[klow];
643 if ( mDxCorr(ilow, jlow) > -1000. ) continue;
644 TMatrixF &mDyCorr = *fLookUpDyCorr[klow];
645 TMatrixF &mDzCorr = *fLookUpDzCorr[klow];
646
647 xref[0]= rRef * TMath::Cos(phiRef);
648 xref[1]= rRef * TMath::Sin(phiRef);
649 xref[2]= zRef;
650
651 FindClosestPosition(ir,iz,iPhi, xref, x2);
652
653 GetDistortion(x2,roc,dx);
654
655 mDxCorr(ilow, jlow) = -dx[0];
656 mDyCorr(ilow, jlow) = -dx[1];
657 mDzCorr(ilow, jlow) = -dx[2];
658
e802250c 659// printf("%3d %3d %3d\n",iPhi, ir, iz);
660// printf("%3d %3d %3d\n",klow, ilow, jlow);
661// printf("x2: %.5f %.5f %.5f\n", x2[0], x2[1], x2[2]);
662// printf("x2d: %.5f %.5f %.5f\n", x2[0]+dx[0], x2[1]+dx[1], x2[2]+dx[2]);
663// printf("xref: %.5f %.5f %.5f\n", xref[0], xref[1], xref[2]);
664// printf("xrd: %.5f %.5f %.5f\n", x2[0]+dx[0]-xref[0], x2[1]+dx[1]-xref[1], x2[2]+dx[2]-xref[2]);
665// printf("phid: %.5f %.5f %.5f\n", phid,rd,zd);
666// printf("phir: %.5f %.5f %.5f\n", phiRef,rRef,zRef);
667// printf("\n");
fd245273 668 }
669 }
670 }
671
672 // fill remaining empty bins
673 // The last ein first phi bin entries must be identical, fill those first
674 {
675 TMatrixF &mDxCorr = *fLookUpDxCorr[0];
676 TMatrixF &mDyCorr = *fLookUpDyCorr[0];
677 TMatrixF &mDzCorr = *fLookUpDzCorr[0];
678
679 TMatrixF &mDxCorr2 = *fLookUpDxCorr[fNPhi-1];
680 TMatrixF &mDyCorr2 = *fLookUpDyCorr[fNPhi-1];
681 TMatrixF &mDzCorr2 = *fLookUpDzCorr[fNPhi-1];
682
683 for (Int_t ir=0; ir<fNR; ++ir){
684 for (Int_t iz=0; iz<fNZ; ++iz){
685 mDxCorr2(ir,iz) = mDxCorr(ir,iz);
686 mDyCorr2(ir,iz) = mDyCorr(ir,iz);
687 mDzCorr2(ir,iz) = mDzCorr(ir,iz);
688 }
689 }
690 }
691
692 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
693 TMatrixF &mDxCorr = *fLookUpDxCorr[iPhi];
694 TMatrixF &mDyCorr = *fLookUpDyCorr[iPhi];
695 TMatrixF &mDzCorr = *fLookUpDzCorr[iPhi];
696
697 Double_t phi=fLimitsPhi(iPhi);
698 for (Int_t ir=0; ir<fNR; ++ir){
699 Double_t r=fLimitsR(ir);
700 x[0]=r * TMath::Cos(phi);
701 x[1]=r * TMath::Sin(phi);
702
703 for (Int_t iz=0; iz<fNZ; ++iz){
704 if (mDxCorr(ir,iz) > -999.) continue;
705
706 Double_t z=fLimitsZ(iz);
707 x[2]=z;
708
709 //TODO: change hardcoded value for r>133.?
710 Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
711 if (r>133.) roc+=36;
712 if (z<0) roc+=18;
713
714 // get last point
715 dx[0] = mDxCorr(ir,iz-1);
716 dx[1] = mDyCorr(ir,iz-1);
717 dx[2] = mDzCorr(ir,iz-1);
718
719 xd[0] = x[0]+dx[0];
720 xd[1] = x[1]+dx[1];
721 xd[2] = x[2]+dx[2];
722
723 // get distorted point
724 const Double_t phid = TVector2::Phi_0_2pi(TMath::ATan2(xd[1],xd[0]));
725 const Double_t rd = TMath::Sqrt(xd[0]*xd[0] + xd[1]*xd[1]);
726 const Double_t zd = xd[2];
727
728 Int_t ilow = 0, jlow = 0, klow = 0 ;
729
730 Search( fLimitsR.GetNrows(), fLimitsR.GetMatrixArray(), rd, ilow ) ;
731 Search( fLimitsZ.GetNrows(), fLimitsZ.GetMatrixArray(), zd, jlow ) ;
732 Search( fLimitsPhi.GetNrows(), fLimitsPhi.GetMatrixArray(), phid, klow ) ;
733
734 if ( ilow < 0 ) ilow = 0 ; // check if out of range
735 if ( jlow < 0 ) jlow = 0 ;
736 if ( klow < 0 ) klow = 0 ;
737 if ( ilow >= fLimitsR.GetNrows()) ilow = fLimitsR.GetNrows() - 1;
738 if ( jlow >= fLimitsZ.GetNrows()) jlow = fLimitsZ.GetNrows() - 1;
739 if ( klow >= fLimitsPhi.GetNrows()) klow = fLimitsPhi.GetNrows() - 1;
740
741 FindClosestPosition(ilow,jlow,klow, x, x2);
742
743 GetDistortion(x2,roc,dx);
744
745 mDxCorr(ir, iz) = -dx[0];
746 mDyCorr(ir, iz) = -dx[1];
747 mDzCorr(ir, iz) = -dx[2];
748 }
749 }
750 }
751
752}
753
754//_________________________________________________________________________________________
755void AliTPCCorrectionLookupTable::FindClosestPosition(const Int_t binR, const Int_t binZ, const Int_t binPhi,
756 const Float_t xref[3], Float_t xret[3])
757{
758 //
759 //
760 //
761
762// static TLinearFitter fitx(2,"pol2");
763// static TLinearFitter fity(2,"pol2");
764// static TLinearFitter fitz(2,"pol2");
765 static TLinearFitter fitx(4,"hyp3");
766 static TLinearFitter fity(4,"hyp3");
767 static TLinearFitter fitz(4,"hyp3");
768 fitx.ClearPoints();
769 fity.ClearPoints();
770 fitz.ClearPoints();
771
772 const Int_t nPoints=3;
773 Int_t counter=0;
774 Int_t rMin=binR;
775 Int_t zMin=binZ;
776 Int_t phiMin=binPhi;
777
778 counter=nPoints/2;
779 while (rMin>0 && counter--) --rMin;
780 counter=nPoints/2;
781 while (zMin>0 && counter--) --zMin;
782 counter=nPoints/2;
783 while (phiMin>0 && counter--) --phiMin;
784
785 Int_t rMax = rMin +nPoints;
786 Int_t zMax = zMin +nPoints;
787 Int_t phiMax = phiMin+nPoints;
788
789 while (rMax>=fNR) {--rMin; --rMax;}
790 while (zMax>=fNZ) {--zMin; --zMax;}
791 while (phiMax>=fNPhi) {--phiMin; --phiMax;}
792
793 Float_t x[3] = {0.,0.,0.};
794 Double_t xd[3] = {0.,0.,0.};
795 Float_t dx[3] = {0.,0.,0.};
796
797 for (Int_t iPhi=phiMin; iPhi<phiMax; ++iPhi) {
798 TMatrixF &mDxDist = *fLookUpDxDist[iPhi];
799 TMatrixF &mDyDist = *fLookUpDyDist[iPhi];
800 TMatrixF &mDzDist = *fLookUpDzDist[iPhi];
801
802 Double_t phi=fLimitsPhi(iPhi);
803 for (Int_t ir=rMin; ir<rMax; ++ir){
804 Double_t r=fLimitsR(ir);
805 x[0]=r * TMath::Cos(phi);
806 x[1]=r * TMath::Sin(phi);
807
808 for (Int_t iz=zMin; iz<zMax; ++iz){
809 Double_t z=fLimitsZ(iz);
810 x[2]=z;
811
812 dx[0] = mDxDist(ir,iz);
813 dx[1] = mDyDist(ir,iz);
814 dx[2] = mDzDist(ir,iz);
815
816 xd[0] = x[0]+dx[0];
817 xd[1] = x[1]+dx[1];
818 xd[2] = x[2]+dx[2];
819
820 fitx.AddPoint(xd, x[0]);
821 fity.AddPoint(xd, x[1]);
822 fitz.AddPoint(xd, x[2]);
823 }
824 }
825 }
826
827 fitx.Eval();
828 fity.Eval();
829 fitz.Eval();
830 xret[0] = fitx.GetParameter(0) + fitx.GetParameter(1)*xref[0]
831 + fitx.GetParameter(2)*xref[1]
832 + fitx.GetParameter(3)*xref[2];
833 xret[1] = fity.GetParameter(0) + fity.GetParameter(1)*xref[0]
834 + fity.GetParameter(2)*xref[1]
835 + fity.GetParameter(3)*xref[2];
836 xret[2] = fitz.GetParameter(0) + fitz.GetParameter(1)*xref[0]
837 + fitz.GetParameter(2)*xref[1]
838 + fitz.GetParameter(3)*xref[2];
839// xret[0] = fitx.GetParameter(0) + fitx.GetParameter(1)*xref[0] + fitx.GetParameter(2)*xref[0]*xref[0];
840// xret[1] = fity.GetParameter(0) + fity.GetParameter(1)*xref[1] + fity.GetParameter(2)*xref[1]*xref[1];
841// xret[2] = fitz.GetParameter(0) + fitz.GetParameter(1)*xref[2] + fitz.GetParameter(2)*xref[2]*xref[2];
842}
843