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 #include <TStopwatch.h>
21 #include <TObjArray.h>
23 #include <THashList.h>
25 #include <TLinearFitter.h>
28 #include <AliTPCROC.h>
30 #include "AliTPCCorrection.h"
32 #include "AliTPCCorrectionLookupTable.h"
34 ClassImp(AliTPCCorrectionLookupTable)
36 //_________________________________________________________________________________________
37 AliTPCCorrectionLookupTable::AliTPCCorrectionLookupTable()
42 , fCorrScaleFactor(-1)
43 , fFillCorrection(kTRUE)
58 //_________________________________________________________________________________________
59 AliTPCCorrectionLookupTable::~AliTPCCorrectionLookupTable()
68 //_________________________________________________________________________________________
69 void AliTPCCorrectionLookupTable::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
71 // Get interpolated Distortion
74 GetInterpolation(x,roc,dx,fLookUpDxDist,fLookUpDyDist,fLookUpDzDist);
77 //_________________________________________________________________________________________
78 void AliTPCCorrectionLookupTable::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
80 // Get interplolated correction
82 GetInterpolation(x,roc,dx,fLookUpDxCorr,fLookUpDyCorr,fLookUpDzCorr);
84 if (fCorrScaleFactor>0) {
85 dx[0]*=fCorrScaleFactor;
86 dx[1]*=fCorrScaleFactor;
87 dx[2]*=fCorrScaleFactor;
91 //_________________________________________________________________________________________
92 void AliTPCCorrectionLookupTable::GetInterpolation(const Float_t x[],const Short_t roc,Float_t dx[],
93 TMatrixF **mDx, TMatrixF **mDy, TMatrixF **mDz)
96 // Calculates the correction/distotring from a lookup table
97 // type: 0 = correction
101 // Float_t typeSign=-1;
102 // if (type==1) typeSign=1;
104 Int_t order = 1 ; // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
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]
114 if ( (roc%36) < 18 ) {
115 sign = 1; // (TPC A side)
117 sign = -1; // (TPC C side)
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
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!");
127 // Get the Er and Ephi field integrals plus the integral over Z
128 dx[0] = Interpolate3DTable(order, r, z, phi,
130 fLimitsR.GetMatrixArray(),
131 fLimitsZ.GetMatrixArray(),
132 fLimitsPhi.GetMatrixArray(),
134 dx[1] = Interpolate3DTable(order, r, z, phi,
136 fLimitsR.GetMatrixArray(),
137 fLimitsZ.GetMatrixArray(),
138 fLimitsPhi.GetMatrixArray(),
140 dx[2] = Interpolate3DTable(order, r, z, phi,
142 fLimitsR.GetMatrixArray(),
143 fLimitsZ.GetMatrixArray(),
144 fLimitsPhi.GetMatrixArray(),
148 //_________________________________________________________________________________________
149 void AliTPCCorrectionLookupTable::CreateLookupTable(AliTPCCorrection &tpcCorr, Float_t stepSize/*=5.*/)
152 // create lookup table for all phi,r,z bins
156 AliError("Limits are not set yet. Please use one of the Set..Limits functions first");
165 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
166 CreateLookupTablePhiBin(tpcCorr,iPhi,stepSize);
170 AliInfo(Form("Required time for lookup table creation: %.2f (%.2f) sec. real (cpu)",s.RealTime(),s.CpuTime()));
173 //_________________________________________________________________________________________
174 void AliTPCCorrectionLookupTable::CreateLookupTableSinglePhi(AliTPCCorrection &tpcCorr, Int_t iPhi, Float_t stepSize)
177 // Lookup table for only one phi bin. Can be used for parallel processing
181 AliError("Limits are not set yet. Please use one of the Set..Limits functions first");
189 InitTablesPhiBin(iPhi);
190 CreateLookupTablePhiBin(tpcCorr,iPhi,stepSize);
193 AliInfo(Form("Required time for lookup table creation: %.2f (%.2f) sec. real (cpu)",s.RealTime(),s.CpuTime()));
196 //_________________________________________________________________________________________
197 void AliTPCCorrectionLookupTable::CreateLookupTablePhiBin(AliTPCCorrection &tpcCorr, Int_t iPhi, Float_t stepSize)
203 if (iPhi<0||iPhi>=fNPhi) return;
205 const Float_t delta=stepSize; // 5cm
206 Float_t x[3]={0.,0.,0.};
207 Float_t dx[3]={0.,0.,0.};
209 Double_t phi=fLimitsPhi(iPhi);
211 TMatrixF &mDxDist = *fLookUpDxDist[iPhi];
212 TMatrixF &mDyDist = *fLookUpDyDist[iPhi];
213 TMatrixF &mDzDist = *fLookUpDzDist[iPhi];
215 TMatrixF &mDxCorr = *fLookUpDxCorr[iPhi];
216 TMatrixF &mDyCorr = *fLookUpDyCorr[iPhi];
217 TMatrixF &mDzCorr = *fLookUpDzCorr[iPhi];
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);
224 for (Int_t iz=0; iz<fNZ; ++iz){
225 Double_t z=fLimitsZ(iz);
227 //TODO: change hardcoded value for r>133.?
228 Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
233 tpcCorr.GetDistortionIntegralDz(x,roc,dx,delta);
235 tpcCorr.GetDistortion(x,roc,dx);
236 mDxDist(ir,iz)=dx[0];
237 mDyDist(ir,iz)=dx[1];
238 mDzDist(ir,iz)=dx[2];
240 if (fFillCorrection) {
242 tpcCorr.GetCorrectionIntegralDz(x,roc,dx,delta);
244 tpcCorr.GetCorrection(x,roc,dx);
245 mDxCorr(ir,iz)=dx[0];
246 mDyCorr(ir,iz)=dx[1];
247 mDzCorr(ir,iz)=dx[2];
254 //_________________________________________________________________________________________
255 void AliTPCCorrectionLookupTable::InitTables()
262 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
263 InitTablesPhiBin(iPhi);
267 //_________________________________________________________________________________________
268 void AliTPCCorrectionLookupTable::CreateLookupTableFromResidualDistortion(THn &resDist)
271 // create lookup table from residual distortions stored in a 3d histogram
272 // assume dimensions are r, phi, z
275 AliError("Limits are not set yet. Please use one of the Set..Limits functions first");
282 Double_t x[3]={0.,0.,0.};
284 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
285 const Double_t phi=fLimitsPhi(iPhi);
288 TMatrixF &mDxDist = *fLookUpDxDist[iPhi];
289 TMatrixF &mDyDist = *fLookUpDyDist[iPhi];
290 TMatrixF &mDzDist = *fLookUpDzDist[iPhi];
292 TMatrixF &mDxCorr = *fLookUpDxCorr[iPhi];
293 TMatrixF &mDyCorr = *fLookUpDyCorr[iPhi];
294 TMatrixF &mDzCorr = *fLookUpDzCorr[iPhi];
296 for (Int_t ir=0; ir<fNR; ++ir){
297 const Double_t r=fLimitsR(ir);
300 for (Int_t iz=0; iz<fNZ; ++iz){
301 const Double_t z=fLimitsZ(iz);
304 const Double_t drphi = resDist.GetBinContent(resDist.GetBin(x));
305 Double_t dx[3]={0.,drphi,0.};
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;
313 mDxDist(ir,iz)=dx[0];
314 mDyDist(ir,iz)=dx[1];
315 mDzDist(ir,iz)=dx[2];
317 mDxCorr(ir,iz)=-dx[0];
318 mDyCorr(ir,iz)=-dx[1];
319 mDzCorr(ir,iz)=-dx[2];
325 //_________________________________________________________________________________________
326 void AliTPCCorrectionLookupTable::CreateResidual(AliTPCCorrection *distortion, AliTPCCorrection* correction)
329 // create lookup table from residual distortions calculated from distorted - correction
335 Float_t x[3]={0.,0.,0.};
337 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
338 const Double_t phi=fLimitsPhi(iPhi);
340 TMatrixF &mDxDist = *fLookUpDxDist[iPhi];
341 TMatrixF &mDyDist = *fLookUpDyDist[iPhi];
342 TMatrixF &mDzDist = *fLookUpDzDist[iPhi];
344 TMatrixF &mDxCorr = *fLookUpDxCorr[iPhi];
345 TMatrixF &mDyCorr = *fLookUpDyCorr[iPhi];
346 TMatrixF &mDzCorr = *fLookUpDzCorr[iPhi];
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);
353 for (Int_t iz=0; iz<fNZ; ++iz){
354 const Double_t z=fLimitsZ(iz);
358 Float_t xdc[3]={x[0], x[1], x[2]};
360 Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
364 //get residual distortion
365 distortion->DistortPoint(xdc, roc);
366 correction->CorrectPoint(xdc, roc);
367 Float_t dx[3]={xdc[0]-x[0], xdc[1]-x[1], xdc[2]-x[2]};
369 mDxDist(ir,iz)=dx[0];
370 mDyDist(ir,iz)=dx[1];
371 mDzDist(ir,iz)=dx[2];
373 mDxCorr(ir,iz)=-dx[0];
374 mDyCorr(ir,iz)=-dx[1];
375 mDzCorr(ir,iz)=-dx[2];
381 //_________________________________________________________________________________________
382 void AliTPCCorrectionLookupTable::InitTablesPhiBin(Int_t iPhi)
388 // check if already initialised
389 if (iPhi<0||iPhi>=fNPhi) return;
390 if (fLookUpDxCorr[iPhi]) return;
392 fLookUpDxCorr[iPhi] = new TMatrixF(fNR,fNZ);
393 fLookUpDyCorr[iPhi] = new TMatrixF(fNR,fNZ);
394 fLookUpDzCorr[iPhi] = new TMatrixF(fNR,fNZ);
396 fLookUpDxDist[iPhi] = new TMatrixF(fNR,fNZ);
397 fLookUpDyDist[iPhi] = new TMatrixF(fNR,fNZ);
398 fLookUpDzDist[iPhi] = new TMatrixF(fNR,fNZ);
401 //_________________________________________________________________________________________
402 void AliTPCCorrectionLookupTable::InitTableArrays()
408 // needs to be before the table creation to set the limits
409 SetupDefaultLimits();
411 fLookUpDxCorr = new TMatrixF*[fNPhi];
412 fLookUpDyCorr = new TMatrixF*[fNPhi];
413 fLookUpDzCorr = new TMatrixF*[fNPhi];
415 fLookUpDxDist = new TMatrixF*[fNPhi];
416 fLookUpDyDist = new TMatrixF*[fNPhi];
417 fLookUpDzDist = new TMatrixF*[fNPhi];
419 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
420 fLookUpDxCorr[iPhi] = 0x0;
421 fLookUpDyCorr[iPhi] = 0x0;
422 fLookUpDzCorr[iPhi] = 0x0;
424 fLookUpDxDist[iPhi] = 0x0;
425 fLookUpDyDist[iPhi] = 0x0;
426 fLookUpDzDist[iPhi] = 0x0;
430 //_________________________________________________________________________________________
431 void AliTPCCorrectionLookupTable::ResetTables()
434 // Reset the lookup tables
437 if (!fLookUpDxCorr) return;
439 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
440 delete fLookUpDxCorr[iPhi];
441 delete fLookUpDyCorr[iPhi];
442 delete fLookUpDzCorr[iPhi];
444 delete fLookUpDxDist[iPhi];
445 delete fLookUpDyDist[iPhi];
446 delete fLookUpDzDist[iPhi];
449 delete [] fLookUpDxCorr;
450 delete [] fLookUpDyCorr;
451 delete [] fLookUpDzCorr;
453 delete [] fLookUpDxDist;
454 delete [] fLookUpDyDist;
455 delete [] fLookUpDzDist;
466 //_________________________________________________________________________________________
467 void AliTPCCorrectionLookupTable::SetupDefaultLimits()
470 // Set default limits for tables
476 fLimitsR. ResizeTo(fNR);
477 fLimitsPhi.ResizeTo(fNPhi);
478 fLimitsZ. ResizeTo(fNZ);
479 fLimitsR. SetElements(fgkRList);
480 fLimitsPhi.SetElements(fgkPhiList);
481 fLimitsZ. SetElements(fgkZList);
484 //_________________________________________________________________________________________
485 void AliTPCCorrectionLookupTable::ResetLimits()
491 fLimitsR. ResizeTo(1);
492 fLimitsPhi.ResizeTo(1);
493 fLimitsZ. ResizeTo(1);
496 //_________________________________________________________________________________________
497 void AliTPCCorrectionLookupTable::MergePhiTables(const char* files)
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
505 ResetLimits(); // use limits from the first file assuming they are all the same
507 TString sfiles=gSystem->GetFromPipe(Form("ls %s",files));
508 TObjArray *arrFiles=sfiles.Tokenize("\n");
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()));
515 AliError(Form("first object in file '%s' is not of type AliTPCCorrectionLookupTable!",f.GetName()));
523 fLimitsR = lt->fLimitsR;
524 fLimitsZ = lt->fLimitsZ;
525 fLimitsPhi = lt->fLimitsPhi;
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()));
534 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi) {
535 if (!lt->fLookUpDxCorr[iPhi]) continue;
537 AliInfo(Form("Adding phi bin '%d' from file '%s'",iPhi,f.GetName()));
539 InitTablesPhiBin(iPhi);
541 *fLookUpDxDist[iPhi] = *lt->fLookUpDxDist[iPhi];
542 *fLookUpDyDist[iPhi] = *lt->fLookUpDyDist[iPhi];
543 *fLookUpDzDist[iPhi] = *lt->fLookUpDzDist[iPhi];
545 *fLookUpDxCorr[iPhi] = *lt->fLookUpDxCorr[iPhi];
546 *fLookUpDyCorr[iPhi] = *lt->fLookUpDyCorr[iPhi];
547 *fLookUpDzCorr[iPhi] = *lt->fLookUpDzCorr[iPhi];
552 //check of all phi bins are initialised
553 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi) {
554 if (!fLookUpDxCorr[iPhi]) {
555 AliFatal(Form("Phi bin '%d' not initialised from files!",iPhi));
562 //_________________________________________________________________________________________
563 void AliTPCCorrectionLookupTable::BuildExactInverse()
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
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.};
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];
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.;
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];
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);
604 for (Int_t iz=0; iz<fNZ; ++iz){
605 Double_t z=fLimitsZ(iz);
608 //TODO: change hardcoded value for r>133.?
609 Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
613 dx[0] = mDxDist(ir,iz);
614 dx[1] = mDyDist(ir,iz);
615 dx[2] = mDzDist(ir,iz);
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];
625 Int_t ilow = 0, jlow = 0, klow = 0 ;
627 Search( fLimitsR.GetNrows(), fLimitsR.GetMatrixArray(), rd, ilow ) ;
628 Search( fLimitsZ.GetNrows(), fLimitsZ.GetMatrixArray(), zd, jlow ) ;
629 Search( fLimitsPhi.GetNrows(), fLimitsPhi.GetMatrixArray(), phid, klow ) ;
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;
638 const Double_t phiRef = fLimitsPhi[klow];
639 const Double_t rRef = fLimitsR[ilow];
640 const Double_t zRef = fLimitsZ[jlow];
642 TMatrixF &mDxCorr = *fLookUpDxCorr[klow];
643 if ( mDxCorr(ilow, jlow) > -1000. ) continue;
644 TMatrixF &mDyCorr = *fLookUpDyCorr[klow];
645 TMatrixF &mDzCorr = *fLookUpDzCorr[klow];
647 xref[0]= rRef * TMath::Cos(phiRef);
648 xref[1]= rRef * TMath::Sin(phiRef);
651 FindClosestPosition(ir,iz,iPhi, xref, x2);
653 GetDistortion(x2,roc,dx);
655 mDxCorr(ilow, jlow) = -dx[0];
656 mDyCorr(ilow, jlow) = -dx[1];
657 mDzCorr(ilow, jlow) = -dx[2];
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);
672 // fill remaining empty bins
673 // The last ein first phi bin entries must be identical, fill those first
675 TMatrixF &mDxCorr = *fLookUpDxCorr[0];
676 TMatrixF &mDyCorr = *fLookUpDyCorr[0];
677 TMatrixF &mDzCorr = *fLookUpDzCorr[0];
679 TMatrixF &mDxCorr2 = *fLookUpDxCorr[fNPhi-1];
680 TMatrixF &mDyCorr2 = *fLookUpDyCorr[fNPhi-1];
681 TMatrixF &mDzCorr2 = *fLookUpDzCorr[fNPhi-1];
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);
692 for (Int_t iPhi=0; iPhi<fNPhi; ++iPhi){
693 TMatrixF &mDxCorr = *fLookUpDxCorr[iPhi];
694 TMatrixF &mDyCorr = *fLookUpDyCorr[iPhi];
695 TMatrixF &mDzCorr = *fLookUpDzCorr[iPhi];
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);
703 for (Int_t iz=0; iz<fNZ; ++iz){
704 if (mDxCorr(ir,iz) > -999.) continue;
706 Double_t z=fLimitsZ(iz);
709 //TODO: change hardcoded value for r>133.?
710 Int_t roc=TMath::Nint(phi*TMath::RadToDeg()/20.)%18;
715 dx[0] = mDxCorr(ir,iz-1);
716 dx[1] = mDyCorr(ir,iz-1);
717 dx[2] = mDzCorr(ir,iz-1);
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];
728 Int_t ilow = 0, jlow = 0, klow = 0 ;
730 Search( fLimitsR.GetNrows(), fLimitsR.GetMatrixArray(), rd, ilow ) ;
731 Search( fLimitsZ.GetNrows(), fLimitsZ.GetMatrixArray(), zd, jlow ) ;
732 Search( fLimitsPhi.GetNrows(), fLimitsPhi.GetMatrixArray(), phid, klow ) ;
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;
741 FindClosestPosition(ilow,jlow,klow, x, x2);
743 GetDistortion(x2,roc,dx);
745 mDxCorr(ir, iz) = -dx[0];
746 mDyCorr(ir, iz) = -dx[1];
747 mDzCorr(ir, iz) = -dx[2];
754 //_________________________________________________________________________________________
755 void AliTPCCorrectionLookupTable::FindClosestPosition(const Int_t binR, const Int_t binZ, const Int_t binPhi,
756 const Float_t xref[3], Float_t xret[3])
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");
772 const Int_t nPoints=3;
779 while (rMin>0 && counter--) --rMin;
781 while (zMin>0 && counter--) --zMin;
783 while (phiMin>0 && counter--) --phiMin;
785 Int_t rMax = rMin +nPoints;
786 Int_t zMax = zMin +nPoints;
787 Int_t phiMax = phiMin+nPoints;
789 while (rMax>=fNR) {--rMin; --rMax;}
790 while (zMax>=fNZ) {--zMin; --zMax;}
791 while (phiMax>=fNPhi) {--phiMin; --phiMax;}
793 Float_t x[3] = {0.,0.,0.};
794 Double_t xd[3] = {0.,0.,0.};
795 Float_t dx[3] = {0.,0.,0.};
797 for (Int_t iPhi=phiMin; iPhi<phiMax; ++iPhi) {
798 TMatrixF &mDxDist = *fLookUpDxDist[iPhi];
799 TMatrixF &mDyDist = *fLookUpDyDist[iPhi];
800 TMatrixF &mDzDist = *fLookUpDzDist[iPhi];
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);
808 for (Int_t iz=zMin; iz<zMax; ++iz){
809 Double_t z=fLimitsZ(iz);
812 dx[0] = mDxDist(ir,iz);
813 dx[1] = mDyDist(ir,iz);
814 dx[2] = mDzDist(ir,iz);
820 fitx.AddPoint(xd, x[0]);
821 fity.AddPoint(xd, x[1]);
822 fitz.AddPoint(xd, x[2]);
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];