Adding option to skip leading tracks in randomized events
[u/mrichter/AliRoot.git] / TPC / AliTPCCorrection.cxx
CommitLineData
0116859c 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// AliTPCCorrection class //
19// //
20// This class provides a general framework to deal with space point //
21// distortions. An correction class which inherits from here is for example //
22// AliTPCExBBShape or AliTPCExBTwist //
23// //
24// General functions are (for example): //
25// CorrectPoint(x,roc) where x is the vector of inital positions in //
26// cartesian coordinates and roc represents the Read Out chamber number //
27// according to the offline naming convention. The vector x is overwritten //
28// with the corrected coordinates. //
29// //
30// An alternative usage would be CorrectPoint(x,roc,dx), which leaves the //
31// vector x untouched, put returns the distortions via the vector dx //
32// //
33// The class allows "effective Omega Tau" corrections to be shifted to the //
34// single distortion classes. //
35// //
36// Note: This class is normally used via the class AliTPCComposedCorrection //
37// //
38// date: 27/04/2010 //
39// Authors: Magnus Mager, Stefan Rossegger, Jim Thomas //
40////////////////////////////////////////////////////////////////////////////////
be67055b 41#include "Riostream.h"
0116859c 42
43#include <TH2F.h>
44#include <TMath.h>
45#include <TROOT.h>
cf5b0aa0 46#include <TTreeStream.h>
ffab0c37 47#include <TTree.h>
48#include <TFile.h>
e527a1b9 49#include <TTimeStamp.h>
ffab0c37 50#include <AliCDBStorage.h>
51#include <AliCDBId.h>
52#include <AliCDBMetaData.h>
7f4cb119 53#include "TVectorD.h"
54
1b923461 55
56#include "TRandom.h"
57#include "AliExternalTrackParam.h"
58#include "AliTrackPointArray.h"
59#include "TDatabasePDG.h"
60#include "AliTrackerBase.h"
61#include "AliTPCROC.h"
62#include "THnSparse.h"
63
7f4cb119 64#include "TRandom.h"
65#include "AliTPCTransform.h"
66#include "AliTPCcalibDB.h"
67#include "AliTPCExB.h"
68#include "AliTPCCorrection.h"
69#include "AliTPCRecoParam.h"
cf5b0aa0 70
71#include "AliExternalTrackParam.h"
72#include "AliTrackPointArray.h"
73#include "TDatabasePDG.h"
74#include "AliTrackerBase.h"
75#include "AliTPCROC.h"
8b63d99c 76#include "THnSparse.h"
1b923461 77
7f4cb119 78#include "AliTPCLaserTrack.h"
ca58ed4e 79#include "AliESDVertex.h"
80#include "AliVertexerTracks.h"
81#include "TDatabasePDG.h"
82#include "TF1.h"
0116859c 83
84#include "AliTPCCorrection.h"
1b923461 85#include "AliLog.h"
0116859c 86
cf5b0aa0 87ClassImp(AliTPCCorrection)
88
f1817479 89
90TObjArray *AliTPCCorrection::fgVisualCorrection=0;
91// instance of correction for visualization
92
93
0116859c 94// FIXME: the following values should come from the database
b1f0a2a5 95const Double_t AliTPCCorrection::fgkTPCZ0 =249.7; // nominal gating grid position
0116859c 96const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
97const Double_t AliTPCCorrection::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
98const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
99const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
100const Double_t AliTPCCorrection::fgkGG =-70.0; // Gating Grid voltage (volts)
101
102
f1817479 103
0116859c 104// FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
105const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
10684.0, 84.5, 85.0, 85.5, 86.0, 87.0, 88.0,
10790.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
108110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
109130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
110150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
111170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
112190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
113210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
114230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 242.0, 244.0, 246.0, 248.0,
115249.0, 249.5, 250.0, 251.5, 252.0 } ;
116
117const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = {
118-249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
119-240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
120-220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
121-200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
122-180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
123-160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
124-140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
125-120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
126-100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0,
127-80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0,
128-60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0,
129-40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0,
130-20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0,
131-1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0,
132 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
133 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0,
134 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0,
135 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0,
136 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0,
137102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
138122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
139142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
140162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
141182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
142202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0,
143222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0,
144242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ;
145
146
147
148AliTPCCorrection::AliTPCCorrection()
534fd34a 149 : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 150{
151 //
152 // default constructor
153 //
f1817479 154 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
0116859c 155}
156
157AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
534fd34a 158: TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 159{
160 //
161 // default constructor, that set the name and title
162 //
f1817479 163 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
0116859c 164}
165
166AliTPCCorrection::~AliTPCCorrection() {
167 //
168 // virtual destructor
169 //
170}
171
172void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
173 //
174 // Corrects the initial coordinates x (cartesian coordinates)
175 // according to the given effect (inherited classes)
176 // roc represents the TPC read out chamber (offline numbering convention)
177 //
178 Float_t dx[3];
179 GetCorrection(x,roc,dx);
180 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
181}
182
183void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
184 //
185 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
186 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
187 // roc represents the TPC read out chamber (offline numbering convention)
188 //
189 Float_t dx[3];
190 GetCorrection(x,roc,dx);
191 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
192}
193
194void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
195 //
196 // Distorts the initial coordinates x (cartesian coordinates)
197 // according to the given effect (inherited classes)
198 // roc represents the TPC read out chamber (offline numbering convention)
199 //
200 Float_t dx[3];
201 GetDistortion(x,roc,dx);
202 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
203}
204
205void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
206 //
207 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
208 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
209 // roc represents the TPC read out chamber (offline numbering convention)
210 //
211 Float_t dx[3];
212 GetDistortion(x,roc,dx);
213 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
214}
215
216void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
217 //
218 // This function delivers the correction values dx in respect to the inital coordinates x
219 // roc represents the TPC read out chamber (offline numbering convention)
220 // Note: The dx is overwritten by the inherited effectice class ...
221 //
222 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
223}
224
225void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
226 //
227 // This function delivers the distortion values dx in respect to the inital coordinates x
228 // roc represents the TPC read out chamber (offline numbering convention)
229 //
230 GetCorrection(x,roc,dx);
231 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
232}
233
234void AliTPCCorrection::Init() {
235 //
236 // Initialization funtion (not used at the moment)
237 //
238}
239
e527a1b9 240void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
241 //
242 // Update function
243 //
244}
245
0116859c 246void AliTPCCorrection::Print(Option_t* /*option*/) const {
247 //
248 // Print function to check which correction classes are used
249 // option=="d" prints details regarding the setted magnitude
250 // option=="a" prints the C0 and C1 coefficents for calibration purposes
251 //
252 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
253}
254
534fd34a 255void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
0116859c 256 //
257 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
258 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
259 // calibration run
260 //
534fd34a 261 fT1=t1;
262 fT2=t2;
263 //SetOmegaTauT1T2(omegaTau, t1, t2);
0116859c 264}
265
266TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
267 //
268 // Simple plot functionality.
269 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
270 // in respect to position z within the XY plane.
271 // The histogramm has nx times ny entries.
272 //
273
274 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
275 nx,-250.,250.,ny,-250.,250.);
276 Float_t x[3],dx[3];
277 x[2]=z;
278 Int_t roc=z>0.?0:18; // FIXME
279 for (Int_t iy=1;iy<=ny;++iy) {
280 x[1]=h->GetYaxis()->GetBinCenter(iy);
281 for (Int_t ix=1;ix<=nx;++ix) {
282 x[0]=h->GetXaxis()->GetBinCenter(ix);
283 GetCorrection(x,roc,dx);
284 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
285 if (90.<=r0 && r0<=250.) {
286 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
287 h->SetBinContent(ix,iy,r1-r0);
288 }
289 else
290 h->SetBinContent(ix,iy,0.);
291 }
292 }
293 return h;
294}
295
296TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
297 //
298 // Simple plot functionality.
299 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
300 // in respect to position z within the XY plane.
301 // The histogramm has nx times ny entries.
302 //
303
304 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
305 nx,-250.,250.,ny,-250.,250.);
306 Float_t x[3],dx[3];
307 x[2]=z;
308 Int_t roc=z>0.?0:18; // FIXME
309 for (Int_t iy=1;iy<=ny;++iy) {
310 x[1]=h->GetYaxis()->GetBinCenter(iy);
311 for (Int_t ix=1;ix<=nx;++ix) {
312 x[0]=h->GetXaxis()->GetBinCenter(ix);
313 GetCorrection(x,roc,dx);
314 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
315 if (90.<=r0 && r0<=250.) {
316 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
317 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
318
319 Float_t dphi=phi1-phi0;
320 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
321 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
322
323 h->SetBinContent(ix,iy,r0*dphi);
324 }
325 else
326 h->SetBinContent(ix,iy,0.);
327 }
328 }
329 return h;
330}
331
332TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
333 //
334 // Simple plot functionality.
335 // Returns a 2d hisogram which represents the corrections in r direction (dr)
336 // in respect to angle phi within the ZR plane.
337 // The histogramm has nx times ny entries.
338 //
339 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
340 nz,-250.,250.,nr,85.,250.);
341 Float_t x[3],dx[3];
342 for (Int_t ir=1;ir<=nr;++ir) {
343 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
344 x[0]=radius*TMath::Cos(phi);
345 x[1]=radius*TMath::Sin(phi);
346 for (Int_t iz=1;iz<=nz;++iz) {
347 x[2]=h->GetXaxis()->GetBinCenter(iz);
348 Int_t roc=x[2]>0.?0:18; // FIXME
349 GetCorrection(x,roc,dx);
350 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
351 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
352 h->SetBinContent(iz,ir,r1-r0);
353 }
354 }
0116859c 355 return h;
356
357}
358
359TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
360 //
361 // Simple plot functionality.
362 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
363 // in respect to angle phi within the ZR plane.
364 // The histogramm has nx times ny entries.
365 //
366 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
367 nz,-250.,250.,nr,85.,250.);
368 Float_t x[3],dx[3];
369 for (Int_t iz=1;iz<=nz;++iz) {
370 x[2]=h->GetXaxis()->GetBinCenter(iz);
371 Int_t roc=x[2]>0.?0:18; // FIXME
372 for (Int_t ir=1;ir<=nr;++ir) {
373 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
374 x[0]=radius*TMath::Cos(phi);
375 x[1]=radius*TMath::Sin(phi);
376 GetCorrection(x,roc,dx);
377 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
378 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
379 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
380
381 Float_t dphi=phi1-phi0;
382 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
383 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
384
385 h->SetBinContent(iz,ir,r0*dphi);
386 }
387 }
388 return h;
389}
390
391TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
392 const char *xlabel,const char *ylabel,const char *zlabel,
393 Int_t nbinsx,Double_t xlow,Double_t xup,
394 Int_t nbinsy,Double_t ylow,Double_t yup) {
395 //
396 // Helper function to create a 2d histogramm of given size
397 //
398
399 TString hname=name;
400 Int_t i=0;
401 if (gDirectory) {
402 while (gDirectory->FindObject(hname.Data())) {
403 hname =name;
404 hname+="_";
405 hname+=i;
406 ++i;
407 }
408 }
409 TH2F *h=new TH2F(hname.Data(),title,
410 nbinsx,xlow,xup,
411 nbinsy,ylow,yup);
412 h->GetXaxis()->SetTitle(xlabel);
413 h->GetYaxis()->SetTitle(ylabel);
414 h->GetZaxis()->SetTitle(zlabel);
415 h->SetStats(0);
416 return h;
417}
418
419
420// Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
421
422void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
b1f0a2a5 423 const Double_t er[kNZ][kNR], Double_t &erValue ) {
0116859c 424 //
425 // Interpolate table - 2D interpolation
426 //
b1f0a2a5 427 Double_t saveEr[10] ;
0116859c 428
429 Search( kNZ, fgkZList, z, fJLow ) ;
430 Search( kNR, fgkRList, r, fKLow ) ;
431 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
432 if ( fKLow < 0 ) fKLow = 0 ;
433 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
434 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
435
436 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
b1f0a2a5 437 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
0116859c 438 }
b1f0a2a5 439 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
0116859c 440
441}
442
443
444Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
b1f0a2a5 445 const Int_t order, const Double_t x ) {
0116859c 446 //
447 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
448 //
449
450 Double_t y ;
451 if ( order == 2 ) { // Quadratic Interpolation = 2
452 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
453 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
454 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
455 } else { // Linear Interpolation = 1
456 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
457 }
458
459 return (y);
460
461}
462
463
b1f0a2a5 464void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
0116859c 465 //
466 // Search an ordered table by starting at the most recently used point
467 //
468
469 Long_t middle, high ;
470 Int_t ascend = 0, increment = 1 ;
471
472 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
473
474 if ( low < 0 || low > n-1 ) {
475 low = -1 ; high = n ;
476 } else { // Ordered Search phase
477 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
478 if ( low == n-1 ) return ;
479 high = low + 1 ;
480 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
481 low = high ;
482 increment *= 2 ;
483 high = low + increment ;
484 if ( high > n-1 ) { high = n ; break ; }
485 }
486 } else {
487 if ( low == 0 ) { low = -1 ; return ; }
488 high = low - 1 ;
489 while ( (Int_t)( x < xArray[low] ) == ascend ) {
490 high = low ;
491 increment *= 2 ;
492 if ( increment >= high ) { low = -1 ; break ; }
493 else low = high - increment ;
494 }
495 }
496 }
497
498 while ( (high-low) != 1 ) { // Binary Search Phase
499 middle = ( high + low ) / 2 ;
500 if ( (Int_t)( x >= xArray[middle] ) == ascend )
501 low = middle ;
502 else
503 high = middle ;
504 }
505
506 if ( x == xArray[n-1] ) low = n-2 ;
507 if ( x == xArray[0] ) low = 0 ;
508
509}
510
1b923461 511void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, const TMatrixD &chargeDensity,
512 TMatrixD &arrayErOverEz, const Int_t rows,
513 const Int_t columns, const Int_t iterations ) {
514 //
515 // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
516 //
517 // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
518 // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
519 // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
520 // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
521 // you wish to solve Laplaces equation however it should not contain random numbers or you will get
522 // random numbers back as a solution.
523 // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
524 // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
525 // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
526 // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
527 // number of points and solves the problem again. This happens several times until the maximum number
528 // of points has been included in the array.
529 //
530 // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
531 // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
532 //
533 // Original code by Jim Thomas (STAR TPC Collaboration)
534 //
535
536 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
537
538 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
539 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
540 const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
541
542 TMatrixD arrayEr(rows,columns) ;
543 TMatrixD arrayEz(rows,columns) ;
544
545 //Check that number of rows and columns is suitable for a binary expansion
546
547 if ( !IsPowerOfTwo(rows-1) ) {
548 AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
549 return;
550 }
551 if ( !IsPowerOfTwo(columns-1) ) {
552 AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
553 return;
554 }
555
556 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
557 // Allow for different size grid spacing in R and Z directions
558 // Use a binary expansion of the size of the matrix to speed up the solution of the problem
559
560 Int_t iOne = (rows-1)/4 ;
561 Int_t jOne = (columns-1)/4 ;
562 // Solve for N in 2**N, add one.
563 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
564
565 for ( Int_t count = 0 ; count < loops ; count++ ) {
566 // Loop while the matrix expands & the resolution increases.
567
568 Float_t tempGridSizeR = gridSizeR * iOne ;
569 Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
570 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
571
572 // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
573 std::vector<float> coef1(rows) ;
574 std::vector<float> coef2(rows) ;
575
576 for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
577 Float_t radius = fgkIFCRadius + i*gridSizeR ;
578 coef1[i] = 1.0 + tempGridSizeR/(2*radius);
579 coef2[i] = 1.0 - tempGridSizeR/(2*radius);
580 }
581
582 TMatrixD sumChargeDensity(rows,columns) ;
583
584 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
585 Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
586 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
587 if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
588 else {
589 // Add up all enclosed charge density contributions within 1/2 unit in all directions
590 Float_t weight = 0.0 ;
591 Float_t sum = 0.0 ;
592 sumChargeDensity(i,j) = 0.0 ;
593 for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
594 for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
595 if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
596 else
597 weight = 1.0 ;
598 // Note that this is cylindrical geometry
599 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
600 sum += weight*radius ;
601 }
602 }
603 sumChargeDensity(i,j) /= sum ;
604 }
605 sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
606 }
607 }
608
609 for ( Int_t k = 1 ; k <= iterations; k++ ) {
610 // Solve Poisson's Equation
611 // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
612 // as interations increase.
613 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
614 Float_t overRelaxM1 = overRelax - 1.0 ;
615 Float_t overRelaxtempFourth, overRelaxcoef5 ;
616 overRelaxtempFourth = overRelax * tempFourth ;
617 overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
618
619 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
620 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
621
622 arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
623 + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
624 - overRelaxcoef5 * arrayV(i,j)
625 + coef1[i] * arrayV(i+iOne,j)
626 + sumChargeDensity(i,j)
627 ) * overRelaxtempFourth;
628 }
629 }
630
631 if ( k == iterations ) {
632 // After full solution is achieved, copy low resolution solution into higher res array
633 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
634 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
635
636 if ( iOne > 1 ) {
637 arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
638 if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
639 }
640 if ( jOne > 1 ) {
641 arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
642 if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
643 }
644 if ( iOne > 1 && jOne > 1 ) {
645 arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
646 if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
647 if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
648 // Note that this leaves a point at the upper left and lower right corners uninitialized.
649 // -> Not a big deal.
650 }
651
652 }
653 }
654 }
655
656 }
657
658 iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
659 jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
660
661 }
662
663 // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
664 for ( Int_t j = 0 ; j < columns ; j++ ) {
665 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
666 arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
667 arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
668 }
669
670 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
671 for ( Int_t i = 0 ; i < rows ; i++) {
672 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
673 arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
674 arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
675 }
676
677 for ( Int_t i = 0 ; i < rows ; i++) {
678 // Note: go back and compare to old version of this code. See notes below.
679 // JT Test ... attempt to divide by real Ez not Ez to first order
680 for ( Int_t j = 0 ; j < columns ; j++ ) {
681 arrayEz(i,j) += ezField;
682 // This adds back the overall Z gradient of the field (main E field component)
683 }
684 // Warning: (-=) assumes you are using an error potetial without the overall Field included
685 }
686
687 // Integrate Er/Ez from Z to zero
688 for ( Int_t j = 0 ; j < columns ; j++ ) {
689 for ( Int_t i = 0 ; i < rows ; i++ ) {
690 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
691 arrayErOverEz(i,j) = 0.0 ;
692 for ( Int_t k = j ; k < columns ; k++ ) {
693 arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
694 if ( index != 4 ) index = 4; else index = 2 ;
695 }
696 if ( index == 4 ) arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
697 if ( index == 2 ) arrayErOverEz(i,j) +=
698 (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
699 -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
700 if ( j == columns-2 ) arrayErOverEz(i,j) =
701 (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
702 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
703 if ( j == columns-1 ) arrayErOverEz(i,j) = 0.0 ;
704 }
705 }
706
707}
708
709
710
710bda39 711Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1b923461 712 //
713 // Helperfunction: Check if integer is a power of 2
714 //
715 Int_t j = 0;
716 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
717 if ( j == 1 ) return(1) ; // True
718 return(0) ; // False
719}
720
cf5b0aa0 721
b1f0a2a5 722AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
cf5b0aa0 723 //
724 // Fit the track parameters - without and with distortion
725 // 1. Space points in the TPC are simulated along the trajectory
726 // 2. Space points distorted
727 // 3. Fits the non distorted and distroted track to the reference plane at refX
728 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
729 //
730 // trackIn - input track parameters
731 // refX - reference X to fit the track
732 // dir - direction - out=1 or in=-1
733 // pcstream - debug streamer to check the results
734 //
cad404e1 735 // see AliExternalTrackParam.h documentation:
736 // track1.fP[0] - local y (rphi)
737 // track1.fP[1] - z
738 // track1.fP[2] - sinus of local inclination angle
739 // track1.fP[3] - tangent of deep angle
740 // track1.fP[4] - 1/pt
1b923461 741
cf5b0aa0 742 AliTPCROC * roc = AliTPCROC::Instance();
743 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
744 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
745 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
cf5b0aa0 746 const Double_t kMaxSnp = 0.85;
747 const Double_t kSigmaY=0.1;
748 const Double_t kSigmaZ=0.1;
ca58ed4e 749 const Double_t kMaxR=500;
750 const Double_t kMaxZ=500;
cf5b0aa0 751 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
ca58ed4e 752 Int_t npoints1=0;
753 Int_t npoints2=0;
cf5b0aa0 754
be67055b 755 AliExternalTrackParam track(trackIn); //
cf5b0aa0 756 // generate points
757 AliTrackPointArray pointArray0(npoints0);
758 AliTrackPointArray pointArray1(npoints0);
759 Double_t xyz[3];
ca58ed4e 760 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 761 //
762 // simulate the track
763 Int_t npoints=0;
764 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
765 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
ca58ed4e 766 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 767 track.GetXYZ(xyz);
ca58ed4e 768 xyz[0]+=gRandom->Gaus(0,0.00005);
769 xyz[1]+=gRandom->Gaus(0,0.00005);
770 xyz[2]+=gRandom->Gaus(0,0.00005);
771 if (TMath::Abs(track.GetZ())>kMaxZ) break;
772 if (TMath::Abs(track.GetX())>kMaxR) break;
cf5b0aa0 773 AliTrackPoint pIn0; // space point
774 AliTrackPoint pIn1;
ffab0c37 775 Int_t sector= (xyz[2]>0)? 0:18;
cf5b0aa0 776 pointArray0.GetPoint(pIn0,npoints);
777 pointArray1.GetPoint(pIn1,npoints);
778 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
779 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
ffab0c37 780 DistortPoint(distPoint, sector);
cf5b0aa0 781 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
782 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
783 //
784 track.Rotate(alpha);
785 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
786 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
787 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
788 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
789 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
790 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
791 pointArray0.AddPoint(npoints, &pIn0);
792 pointArray1.AddPoint(npoints, &pIn1);
793 npoints++;
794 if (npoints>=npoints0) break;
795 }
7f4cb119 796 if (npoints<npoints0/2) return 0;
cf5b0aa0 797 //
798 // refit track
799 //
800 AliExternalTrackParam *track0=0;
801 AliExternalTrackParam *track1=0;
802 AliTrackPoint point1,point2,point3;
803 if (dir==1) { //make seed inner
804 pointArray0.GetPoint(point1,1);
4486a91f 805 pointArray0.GetPoint(point2,30);
806 pointArray0.GetPoint(point3,60);
cf5b0aa0 807 }
808 if (dir==-1){ //make seed outer
4486a91f 809 pointArray0.GetPoint(point1,npoints-60);
810 pointArray0.GetPoint(point2,npoints-30);
cf5b0aa0 811 pointArray0.GetPoint(point3,npoints-1);
812 }
813 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
814 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
815
cf5b0aa0 816 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
8b63d99c 817 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
cf5b0aa0 818 //
819 AliTrackPoint pIn0;
820 AliTrackPoint pIn1;
821 pointArray0.GetPoint(pIn0,ipoint);
822 pointArray1.GetPoint(pIn1,ipoint);
823 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
824 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
825 //
ca58ed4e 826 if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
827 if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
828 if (TMath::Abs(track0->GetZ())>kMaxZ) break;
829 if (TMath::Abs(track0->GetX())>kMaxR) break;
830 if (TMath::Abs(track1->GetZ())>kMaxZ) break;
831 if (TMath::Abs(track1->GetX())>kMaxR) break;
832
8b63d99c 833 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
cf5b0aa0 834 //
835 Double_t pointPos[2]={0,0};
836 Double_t pointCov[3]={0,0,0};
837 pointPos[0]=prot0.GetY();//local y
838 pointPos[1]=prot0.GetZ();//local z
839 pointCov[0]=prot0.GetCov()[3];//simay^2
840 pointCov[1]=prot0.GetCov()[4];//sigmayz
841 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
ca58ed4e 842 if (!track0->Update(pointPos,pointCov)) break;
cf5b0aa0 843 //
8b63d99c 844 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
845 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
846 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
847
0b736a46 848 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
849 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
cf5b0aa0 850 pointCov[0]=prot1.GetCov()[3];//simay^2
851 pointCov[1]=prot1.GetCov()[4];//sigmayz
852 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
ca58ed4e 853 if (!track1->Update(pointPos,pointCov)) break;
854 npoints1++;
855 npoints2++;
cf5b0aa0 856 }
ca58ed4e 857 if (npoints2<npoints) return 0;
8b63d99c 858 AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
cf5b0aa0 859 track1->Rotate(track0->GetAlpha());
ca58ed4e 860 AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
cf5b0aa0 861
cad404e1 862 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
cf5b0aa0 863 "point0.="<<&pointArray0<< // points
864 "point1.="<<&pointArray1<< // distorted points
865 "trackIn.="<<&track<< // original track
866 "track0.="<<track0<< // fitted track
867 "track1.="<<track1<< // fitted distorted track
868 "\n";
be67055b 869 new(&trackIn) AliExternalTrackParam(*track0);
cf5b0aa0 870 delete track0;
871 return track1;
872}
873
874
ffab0c37 875
876
877
878TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
879 //
880 // create the distortion tree on a mesh with granularity given by step
881 // return the tree with distortions at given position
882 // Map is created on the mesh with given step size
883 //
884 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
885 Float_t xyz[3];
886 for (Double_t x= -250; x<250; x+=step){
887 for (Double_t y= -250; y<250; y+=step){
888 Double_t r = TMath::Sqrt(x*x+y*y);
889 if (r<80) continue;
890 if (r>250) continue;
891 for (Double_t z= -250; z<250; z+=step){
892 Int_t roc=(z>0)?0:18;
893 xyz[0]=x;
894 xyz[1]=y;
895 xyz[2]=z;
896 Double_t phi = TMath::ATan2(y,x);
897 DistortPoint(xyz,roc);
898 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
899 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
900 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
901 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
902 Double_t dx = xyz[0]-x;
903 Double_t dy = xyz[1]-y;
904 Double_t dz = xyz[2]-z;
905 Double_t dr=r1-r;
906 Double_t drphi=(phi1-phi)*r;
907 (*pcstream)<<"distortion"<<
908 "x="<<x<< // original position
909 "y="<<y<<
910 "z="<<z<<
911 "r="<<r<<
912 "phi="<<phi<<
913 "x1="<<xyz[0]<< // distorted position
914 "y1="<<xyz[1]<<
915 "z1="<<xyz[2]<<
916 "r1="<<r1<<
917 "phi1="<<phi1<<
918 //
919 "dx="<<dx<< // delta position
920 "dy="<<dy<<
921 "dz="<<dz<<
922 "dr="<<dr<<
923 "drphi="<<drphi<<
924 "\n";
925 }
926 }
927 }
928 delete pcstream;
929 TFile f(Form("correction%s.root",GetName()));
930 TTree * tree = (TTree*)f.Get("distortion");
931 TTree * tree2= tree->CopyTree("1");
932 tree2->SetName(Form("dist%s",GetName()));
933 tree2->SetDirectory(0);
934 delete tree;
935 return tree2;
936}
937
938
939
be67055b 940
b1f0a2a5 941void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
be67055b 942 //
943 // Make a fit tree:
944 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
945 // calculates partial distortions
946 // Partial distortion is stored in the resulting tree
947 // Output is storred in the file distortion_<dettype>_<partype>.root
948 // Partial distortion is stored with the name given by correction name
949 //
950 //
951 // Parameters of function:
952 // input - input tree
953 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex
954 // ppype - parameter type
955 // corrArray - array with partial corrections
956 // step - skipe entries - if 1 all entries processed - it is slow
957 // debug 0 if debug on also space points dumped - it is slow
b322e06a 958 const Double_t kMaxSnp = 0.85;
959 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
960 // const Double_t kB2C=-0.299792458e-3;
7f4cb119 961 const Int_t kMinEntries=50;
be67055b 962 Double_t phi,theta, snp, mean,rms, entries;
963 tinput->SetBranchAddress("theta",&theta);
964 tinput->SetBranchAddress("phi", &phi);
965 tinput->SetBranchAddress("snp",&snp);
966 tinput->SetBranchAddress("mean",&mean);
967 tinput->SetBranchAddress("rms",&rms);
968 tinput->SetBranchAddress("entries",&entries);
969 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
970 //
971 Int_t nentries=tinput->GetEntries();
972 Int_t ncorr=corrArray->GetEntries();
7f4cb119 973 Double_t corrections[100]={0}; //
be67055b 974 Double_t tPar[5];
975 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
976 Double_t refX=0;
977 Int_t dir=0;
7f4cb119 978 if (dtype==0) {refX=85.; dir=-1;}
979 if (dtype==1) {refX=275.; dir=1;}
980 if (dtype==2) {refX=85.; dir=-1;}
981 if (dtype==3) {refX=360.; dir=-1;}
be67055b 982 //
983 for (Int_t ientry=0; ientry<nentries; ientry+=step){
984 tinput->GetEntry(ientry);
7f4cb119 985 if (TMath::Abs(snp)>kMaxSnp) continue;
be67055b 986 tPar[0]=0;
987 tPar[1]=theta*refX;
988 tPar[2]=snp;
989 tPar[3]=theta;
4486a91f 990 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
8b63d99c 991 Double_t bz=AliTrackerBase::GetBz();
4486a91f 992 if (refX>10. && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
993 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
7f4cb119 994 AliExternalTrackParam track(refX,phi,tPar,cov);
995 Double_t xyz[3];
996 track.GetXYZ(xyz);
997 Int_t id=0;
998 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
0b736a46 999 if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature
be67055b 1000 (*pcstream)<<"fit"<<
8b63d99c 1001 "bz="<<bz<< // magnetic filed used
be67055b 1002 "dtype="<<dtype<< // detector match type
1003 "ptype="<<ptype<< // parameter type
1004 "theta="<<theta<< // theta
1005 "phi="<<phi<< // phi
1006 "snp="<<snp<< // snp
1007 "mean="<<mean<< // mean dist value
1008 "rms="<<rms<< // rms
7f4cb119 1009 "gx="<<xyz[0]<< // global position at reference
1010 "gy="<<xyz[1]<< // global position at reference
1011 "gz="<<xyz[2]<< // global position at reference
1012 "dRrec="<<dRrec<< // delta Radius in reconstruction
1013 "id="<<id<< // track id
be67055b 1014 "entries="<<entries;// number of entries in bin
1015 //
1016 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1017 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1018 corrections[icorr]=0;
1019 if (entries>kMinEntries){
1020 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1021 AliExternalTrackParam *trackOut = 0;
1022 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1023 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
7f4cb119 1024 if (dtype==0) {refX=85.; dir=-1;}
1025 if (dtype==1) {refX=275.; dir=1;}
b1f0a2a5 1026 if (dtype==2) {refX=0; dir=-1;}
7f4cb119 1027 if (dtype==3) {refX=360.; dir=-1;}
b1f0a2a5 1028 //
7f4cb119 1029 if (trackOut){
1030 AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1031 trackOut->Rotate(trackIn.GetAlpha());
1032 trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1033 //
1034 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1035 delete trackOut;
1036 }else{
1037 corrections[icorr]=0;
1038 }
0b736a46 1039 if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature
be67055b 1040 }
7f4cb119 1041 Double_t dRdummy=0;
be67055b 1042 (*pcstream)<<"fit"<<
7f4cb119 1043 Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
1044 Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
1045 // for points it is neccessary
be67055b 1046 }
1047 (*pcstream)<<"fit"<<"\n";
1048 }
1049 delete pcstream;
1050}
1051
1052
1053
7f4cb119 1054void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1055 //
1056 // Make a laser fit tree for global minimization
1057 //
1058 const Double_t cutErrY=0.1;
1059 const Double_t cutErrZ=0.1;
1060 const Double_t kEpsilon=0.00000001;
1061 TVectorD *vecdY=0;
1062 TVectorD *vecdZ=0;
1063 TVectorD *veceY=0;
1064 TVectorD *veceZ=0;
1065 AliTPCLaserTrack *ltr=0;
1066 AliTPCLaserTrack::LoadTracks();
1067 tree->SetBranchAddress("dY.",&vecdY);
1068 tree->SetBranchAddress("dZ.",&vecdZ);
1069 tree->SetBranchAddress("eY.",&veceY);
1070 tree->SetBranchAddress("eZ.",&veceZ);
1071 tree->SetBranchAddress("LTr.",&ltr);
1072 Int_t entries= tree->GetEntries();
1073 TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1074 Double_t bz=AliTrackerBase::GetBz();
1075 //
1076
1077 for (Int_t ientry=0; ientry<entries; ientry++){
1078 tree->GetEntry(ientry);
1079 if (!ltr->GetVecGX()){
1080 ltr->UpdatePoints();
1081 }
1082 TVectorD * delta= (itype==0)? vecdY:vecdZ;
1083 TVectorD * err= (itype==0)? veceY:veceZ;
1084
1085 for (Int_t irow=0; irow<159; irow++){
1086 Int_t nentries = 1000;
1087 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1088 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1089 Int_t dtype=4;
1090 Double_t phi =(*ltr->GetVecPhi())[irow];
1091 Double_t theta =ltr->GetTgl();
1092 Double_t mean=delta->GetMatrixArray()[irow];
1093 Double_t gx=0,gy=0,gz=0;
1094 Double_t snp = (*ltr->GetVecP2())[irow];
1095 Double_t rms = 0.1+err->GetMatrixArray()[irow];
1096 gx = (*ltr->GetVecGX())[irow];
1097 gy = (*ltr->GetVecGY())[irow];
1098 gz = (*ltr->GetVecGZ())[irow];
1099 Int_t bundle= ltr->GetBundle();
1100 Double_t dRrec=0;
1101 //
1102 // get delta R used in reconstruction
1103 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
1104 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1105 const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1106 Double_t xyz0[3]={gx,gy,gz};
1107 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1108 //
1109 // old ExB correction
1110 //
1111 if(recoParam&&recoParam->GetUseExBCorrection()) {
1112 Double_t xyz1[3]={gx,gy,gz};
1113 calib->GetExB()->Correct(xyz0,xyz1);
1114 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1115 dRrec=oldR-newR;
1116 }
1117 if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1118 Float_t xyz1[3]={gx,gy,gz};
1119 Int_t sector=(gz>0)?0:18;
1120 correction->CorrectPoint(xyz1, sector);
1121 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1122 dRrec=oldR-newR;
1123 }
1124
1125
1126 (*pcstream)<<"fit"<<
1127 "bz="<<bz<< // magnetic filed used
1128 "dtype="<<dtype<< // detector match type
1129 "ptype="<<itype<< // parameter type
1130 "theta="<<theta<< // theta
1131 "phi="<<phi<< // phi
1132 "snp="<<snp<< // snp
1133 "mean="<<mean<< // mean dist value
1134 "rms="<<rms<< // rms
1135 "gx="<<gx<< // global position
1136 "gy="<<gy<< // global position
1137 "gz="<<gz<< // global position
1138 "dRrec="<<dRrec<< // delta Radius in reconstruction
1139 "id="<<bundle<< //bundle
1140 "entries="<<nentries;// number of entries in bin
1141 //
1142 //
1143 Double_t ky = TMath::Tan(TMath::ASin(snp));
1144 Int_t ncorr = corrArray->GetEntries();
1145 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
1146 Double_t phi0 = TMath::ATan2(gy,gx);
1147 Double_t distortions[1000]={0};
1148 Double_t distortionsR[1000]={0};
1149 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1150 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1151 Float_t distPoint[3]={gx,gy,gz};
1152 Int_t sector= (gz>0)? 0:18;
1153 if (r0>80){
1154 corr->DistortPoint(distPoint, sector);
1155 }
1b923461 1156 // Double_t value=distPoint[2]-gz;
7f4cb119 1157 if (itype==0){
1158 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1159 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1160 Double_t drphi= r0*(phi1-phi0);
1161 Double_t dr = r1-r0;
1162 distortions[icorr] = drphi-ky*dr;
1163 distortionsR[icorr] = dr;
1164 }
1165 (*pcstream)<<"fit"<<
1166 Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
1167 Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
1168 }
1169 (*pcstream)<<"fit"<<"\n";
1170 }
1171 }
1172 delete pcstream;
1173}
1174
1175
be67055b 1176
b1f0a2a5 1177void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
8b63d99c 1178 //
1179 // make a distortion map out ou fthe residual histogram
1180 // Results are written to the debug streamer - pcstream
1181 // Parameters:
1182 // his0 - input (4D) residual histogram
1183 // pcstream - file to write the tree
1184 // run - run number
1185 // marian.ivanov@cern.ch
1186 const Int_t kMinEntries=50;
1187 Int_t nbins1=his0->GetAxis(1)->GetNbins();
1188 Int_t first1=his0->GetAxis(1)->GetFirst();
1189 Int_t last1 =his0->GetAxis(1)->GetLast();
1190 //
1191 Double_t bz=AliTrackerBase::GetBz();
1192 Int_t idim[4]={0,1,2,3};
1193 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
1194 //
1195 his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1196 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1197 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
1198 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
1199 Int_t first3 = his1->GetAxis(3)->GetFirst();
1200 Int_t last3 = his1->GetAxis(3)->GetLast();
1201 //
1202
1203 for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){ // axis 3 - local angle
1204 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1205 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1206 if (ibin3<first3) {
1207 his1->GetAxis(3)->SetRangeUser(-1,1);
1208 x3=0;
1209 }
1210 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
1211 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1212 Int_t first2 = his3->GetAxis(2)->GetFirst();
1213 Int_t last2 = his3->GetAxis(2)->GetLast();
1214 //
1215 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1216 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1217 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1218 TH1 * hisDelta = his3->Projection(0);
1219 //
1220 Double_t entries = hisDelta->GetEntries();
1221 Double_t mean=0, rms=0;
1222 if (entries>kMinEntries){
1223 mean = hisDelta->GetMean();
1224 rms = hisDelta->GetRMS();
1225 }
1226 (*pcstream)<<hname<<
1227 "run="<<run<<
1228 "bz="<<bz<<
1229 "theta="<<x1<<
1230 "phi="<<x2<<
1231 "snp="<<x3<<
1232 "entries="<<entries<<
1233 "mean="<<mean<<
1234 "rms="<<rms<<
1235 "\n";
1236 delete hisDelta;
1237 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1238 }
1239 delete his3;
1240 }
1241 delete his1;
1242 }
1243}
1244
1245
1246
1247
1248
ffab0c37 1249void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1250 //
1251 // Store object in the OCDB
1252 // By default the object is stored in the current directory
1253 // default comment consit of user name and the date
1254 //
1255 TString ocdbStorage="";
1256 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1257 AliCDBMetaData *metaData= new AliCDBMetaData();
1258 metaData->SetObjectClassName("AliTPCCorrection");
1259 metaData->SetResponsible("Marian Ivanov");
1260 metaData->SetBeamPeriod(1);
1261 metaData->SetAliRootVersion("05-25-01"); //root version
1262 TString userName=gSystem->GetFromPipe("echo $USER");
1263 TString date=gSystem->GetFromPipe("date");
1264
1265 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1266 if (comment) metaData->SetComment(comment);
1267 AliCDBId* id1=NULL;
1268 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1269 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1270 gStorage->Put(this, (*id1), metaData);
1271}
1272
ca58ed4e 1273
7d85e147 1274void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
ca58ed4e 1275
1276 AliVertexerTracks *vertexer = new AliVertexerTracks(5);// 5kGaus
1277
1278 TObjArray ATrk; // Original Track array of Aside
1279 TObjArray dATrk; // Distorted Track array of A side
1280 UShort_t *AId = new UShort_t[nTracks]; // A side Track ID
1281 TObjArray CTrk;
1282 TObjArray dCTrk;
1283 UShort_t *CId = new UShort_t [nTracks];
1284 Int_t ID=0;
1285 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
7d85e147 1286 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
ca58ed4e 1287 fpt.SetParameters(7.24,0.120);
1288 fpt.SetNpx(10000);
1289 for(Int_t nt=0; nt<nTracks; nt++){
1290 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
7d85e147 1291 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
ca58ed4e 1292 Double_t pt = fpt.GetRandom();// momentum for f1
1293 Short_t sign=1;
1294 if(gRandom->Rndm() < 0.5){
1295 sign =1;
1296 }else{
1297 sign=-1;
1298 }
1299
1300 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1301 Double_t pxyz[3];
1302 pxyz[0]=pt*TMath::Cos(phi);
1303 pxyz[1]=pt*TMath::Sin(phi);
1304 pxyz[2]=pt*TMath::Tan(theta);
1305 Double_t cv[21]={0};
1306 AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
1307
1308 Double_t refX=1.;
1309 Int_t dir=-1;
1310 AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir, NULL);
1311 if (!td) continue;
1312 if (pcstream) (*pcstream)<<"track"<<
1313 "eta="<<eta<<
1314 "theta="<<theta<<
1315 "tOrig.="<<t<<
1316 "td.="<<td<<
1317 "\n";
7d85e147 1318 if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
ca58ed4e 1319 if (td){
1320 dATrk.AddLast(td);
1321 ATrk.AddLast(t);
1322 Int_t nn=ATrk.GetEntriesFast();
1323 AId[nn]=ID;
1324 }
7d85e147 1325 }else if(( eta<-0.07 )&&( eta>-etaCuts )){
ca58ed4e 1326 if (td){
1327 dCTrk.AddLast(td);
1328 CTrk.AddLast(t);
1329 Int_t nn=CTrk.GetEntriesFast();
1330 CId[nn]=ID;
1331 }
1332 }
1333 ID++;
1334 }// end of track loop
1335
1336 vertexer->SetTPCMode();
1337 vertexer->SetConstraintOff();
1338
1339 aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dATrk,AId));
1340 avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&ATrk,AId));
1341 cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dCTrk,CId));
1342 cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&CTrk,CId));
1343 if (pcstream) (*pcstream)<<"vertex"<<
1344 "x="<<orgVertex[0]<<
1345 "y="<<orgVertex[1]<<
1346 "z="<<orgVertex[2]<<
1347 "av.="<<&aV<< // distorted vertex A side
1348 "cv.="<<&cV<< // distroted vertex C side
1349 "avO.="<<&avOrg<< // original vertex A side
1350 "cvO.="<<&cvOrg<<
1351 "\n";
1352 delete []AId;
1353 delete []CId;
1354}
f1817479 1355
1356void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
1357 //
1358 // make correction available for visualization using
1359 // TFormula, TFX and TTree::Draw
1360 // important in order to check corrections and also compute dervied variables
1361 // e.g correction partial derivatives
1362 //
1363 // NOTE - class is not owner of correction
1364 //
1365 if (!fgVisualCorrection) fgVisualCorrection=new TObjArray;
1366 fgVisualCorrection->AddAt(corr, position);
1367}
1368
1369
1370
1371Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t localX, Double_t kZ, Int_t axisType, Int_t corrType){
1372 //
1373 // calculate the correction at given position - check the geffCorr
1374 //
1375 if (!fgVisualCorrection) return 0;
1376 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1377 if (!corr) return 0;
1378 Double_t phi=sector*TMath::Pi()/9.;
1379 Double_t gx = localX*TMath::Cos(phi);
1380 Double_t gy = localX*TMath::Sin(phi);
1381 Double_t gz = localX*kZ;
1382 Int_t nsector=(gz>0) ? 0:18;
1383 //
1384 //
1385 //
1386 Float_t distPoint[3]={gx,gy,gz};
1387 corr->DistortPoint(distPoint, nsector);
1388 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1389 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1390 Double_t phi0=TMath::ATan2(gy,gx);
1391 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1392 if (axisType==0) return r1-r0;
1393 if (axisType==1) return (phi1-phi0)*r0;
1394 if (axisType==2) return distPoint[2]-gz;
1395 return phi1-phi0;
1396}
1397
1398Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
1399 //
1400 // return correction at given x,y,z
1401 //
1402 if (!fgVisualCorrection) return 0;
1403 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1404 if (!corr) return 0;
1405 Double_t phi0= TMath::ATan2(gy,gx);
1406 Int_t nsector=(gz>0) ? 0:18;
1407 Float_t distPoint[3]={gx,gy,gz};
1408 corr->DistortPoint(distPoint, nsector);
1409 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1410 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1411 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1412 if (axisType==0) return r1-r0;
1413 if (axisType==1) return (phi1-phi0)*r0;
1414 if (axisType==2) return distPoint[2]-gz;
1415 return phi1-phi0;
1416}