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