]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCCorrection.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
adding the tansformation names to the trees
[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////////////////////////////////////////////////////////////////////////////////
41
42#include <TH2F.h>
43#include <TMath.h>
44#include <TROOT.h>
cf5b0aa0 45#include <TTreeStream.h>
46
47#include "AliExternalTrackParam.h"
48#include "AliTrackPointArray.h"
49#include "TDatabasePDG.h"
50#include "AliTrackerBase.h"
51#include "AliTPCROC.h"
52
0116859c 53
54#include "AliTPCCorrection.h"
55
cf5b0aa0 56ClassImp(AliTPCCorrection)
57
0116859c 58// FIXME: the following values should come from the database
59const Double_t AliTPCCorrection::fgkTPC_Z0 =249.7; // nominal gating grid position
60const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
61const Double_t AliTPCCorrection::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
62const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
63const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
64const Double_t AliTPCCorrection::fgkGG =-70.0; // Gating Grid voltage (volts)
65
66
67// FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
68const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
6984.0, 84.5, 85.0, 85.5, 86.0, 87.0, 88.0,
7090.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
71110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
72130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
73150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
74170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
75190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
76210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
77230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 242.0, 244.0, 246.0, 248.0,
78249.0, 249.5, 250.0, 251.5, 252.0 } ;
79
80const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = {
81-249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
82-240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
83-220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
84-200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
85-180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
86-160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
87-140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
88-120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
89-100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0,
90-80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0,
91-60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0,
92-40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0,
93-20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0,
94-1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0,
95 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
96 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0,
97 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0,
98 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0,
99 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0,
100102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
101122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
102142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
103162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
104182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
105202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0,
106222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0,
107242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ;
108
109
110
111AliTPCCorrection::AliTPCCorrection()
112 : TNamed("correction_unity","unity"),fJLow(0),fKLow(0)
113{
114 //
115 // default constructor
116 //
117}
118
119AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
120 : TNamed(name,title),fJLow(0),fKLow(0)
121{
122 //
123 // default constructor, that set the name and title
124 //
125}
126
127AliTPCCorrection::~AliTPCCorrection() {
128 //
129 // virtual destructor
130 //
131}
132
133void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
134 //
135 // Corrects the initial coordinates x (cartesian coordinates)
136 // according to the given effect (inherited classes)
137 // roc represents the TPC read out chamber (offline numbering convention)
138 //
139 Float_t dx[3];
140 GetCorrection(x,roc,dx);
141 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
142}
143
144void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
145 //
146 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
147 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
148 // roc represents the TPC read out chamber (offline numbering convention)
149 //
150 Float_t dx[3];
151 GetCorrection(x,roc,dx);
152 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
153}
154
155void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
156 //
157 // Distorts the initial coordinates x (cartesian coordinates)
158 // according to the given effect (inherited classes)
159 // roc represents the TPC read out chamber (offline numbering convention)
160 //
161 Float_t dx[3];
162 GetDistortion(x,roc,dx);
163 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
164}
165
166void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
167 //
168 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
169 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
170 // roc represents the TPC read out chamber (offline numbering convention)
171 //
172 Float_t dx[3];
173 GetDistortion(x,roc,dx);
174 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
175}
176
177void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
178 //
179 // This function delivers the correction values dx in respect to the inital coordinates x
180 // roc represents the TPC read out chamber (offline numbering convention)
181 // Note: The dx is overwritten by the inherited effectice class ...
182 //
183 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
184}
185
186void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
187 //
188 // This function delivers the distortion values dx in respect to the inital coordinates x
189 // roc represents the TPC read out chamber (offline numbering convention)
190 //
191 GetCorrection(x,roc,dx);
192 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
193}
194
195void AliTPCCorrection::Init() {
196 //
197 // Initialization funtion (not used at the moment)
198 //
199}
200
201void AliTPCCorrection::Print(Option_t* /*option*/) const {
202 //
203 // Print function to check which correction classes are used
204 // option=="d" prints details regarding the setted magnitude
205 // option=="a" prints the C0 and C1 coefficents for calibration purposes
206 //
207 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
208}
209
210void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t /*t1*/,Float_t /*t2*/) {
211 //
212 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
213 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
214 // calibration run
215 //
216 // SetOmegaTauT1T2(omegaTau, t1, t2);
217}
218
219TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
220 //
221 // Simple plot functionality.
222 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
223 // in respect to position z within the XY plane.
224 // The histogramm has nx times ny entries.
225 //
226
227 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
228 nx,-250.,250.,ny,-250.,250.);
229 Float_t x[3],dx[3];
230 x[2]=z;
231 Int_t roc=z>0.?0:18; // FIXME
232 for (Int_t iy=1;iy<=ny;++iy) {
233 x[1]=h->GetYaxis()->GetBinCenter(iy);
234 for (Int_t ix=1;ix<=nx;++ix) {
235 x[0]=h->GetXaxis()->GetBinCenter(ix);
236 GetCorrection(x,roc,dx);
237 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
238 if (90.<=r0 && r0<=250.) {
239 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
240 h->SetBinContent(ix,iy,r1-r0);
241 }
242 else
243 h->SetBinContent(ix,iy,0.);
244 }
245 }
246 return h;
247}
248
249TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
250 //
251 // Simple plot functionality.
252 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
253 // in respect to position z within the XY plane.
254 // The histogramm has nx times ny entries.
255 //
256
257 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
258 nx,-250.,250.,ny,-250.,250.);
259 Float_t x[3],dx[3];
260 x[2]=z;
261 Int_t roc=z>0.?0:18; // FIXME
262 for (Int_t iy=1;iy<=ny;++iy) {
263 x[1]=h->GetYaxis()->GetBinCenter(iy);
264 for (Int_t ix=1;ix<=nx;++ix) {
265 x[0]=h->GetXaxis()->GetBinCenter(ix);
266 GetCorrection(x,roc,dx);
267 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
268 if (90.<=r0 && r0<=250.) {
269 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
270 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
271
272 Float_t dphi=phi1-phi0;
273 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
274 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
275
276 h->SetBinContent(ix,iy,r0*dphi);
277 }
278 else
279 h->SetBinContent(ix,iy,0.);
280 }
281 }
282 return h;
283}
284
285TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
286 //
287 // Simple plot functionality.
288 // Returns a 2d hisogram which represents the corrections in r direction (dr)
289 // in respect to angle phi within the ZR plane.
290 // The histogramm has nx times ny entries.
291 //
292 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
293 nz,-250.,250.,nr,85.,250.);
294 Float_t x[3],dx[3];
295 for (Int_t ir=1;ir<=nr;++ir) {
296 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
297 x[0]=radius*TMath::Cos(phi);
298 x[1]=radius*TMath::Sin(phi);
299 for (Int_t iz=1;iz<=nz;++iz) {
300 x[2]=h->GetXaxis()->GetBinCenter(iz);
301 Int_t roc=x[2]>0.?0:18; // FIXME
302 GetCorrection(x,roc,dx);
303 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
304 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
305 h->SetBinContent(iz,ir,r1-r0);
306 }
307 }
308 printf("SDF\n");
309 return h;
310
311}
312
313TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
314 //
315 // Simple plot functionality.
316 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
317 // in respect to angle phi within the ZR plane.
318 // The histogramm has nx times ny entries.
319 //
320 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
321 nz,-250.,250.,nr,85.,250.);
322 Float_t x[3],dx[3];
323 for (Int_t iz=1;iz<=nz;++iz) {
324 x[2]=h->GetXaxis()->GetBinCenter(iz);
325 Int_t roc=x[2]>0.?0:18; // FIXME
326 for (Int_t ir=1;ir<=nr;++ir) {
327 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
328 x[0]=radius*TMath::Cos(phi);
329 x[1]=radius*TMath::Sin(phi);
330 GetCorrection(x,roc,dx);
331 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
332 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
333 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
334
335 Float_t dphi=phi1-phi0;
336 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
337 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
338
339 h->SetBinContent(iz,ir,r0*dphi);
340 }
341 }
342 return h;
343}
344
345TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
346 const char *xlabel,const char *ylabel,const char *zlabel,
347 Int_t nbinsx,Double_t xlow,Double_t xup,
348 Int_t nbinsy,Double_t ylow,Double_t yup) {
349 //
350 // Helper function to create a 2d histogramm of given size
351 //
352
353 TString hname=name;
354 Int_t i=0;
355 if (gDirectory) {
356 while (gDirectory->FindObject(hname.Data())) {
357 hname =name;
358 hname+="_";
359 hname+=i;
360 ++i;
361 }
362 }
363 TH2F *h=new TH2F(hname.Data(),title,
364 nbinsx,xlow,xup,
365 nbinsy,ylow,yup);
366 h->GetXaxis()->SetTitle(xlabel);
367 h->GetYaxis()->SetTitle(ylabel);
368 h->GetZaxis()->SetTitle(zlabel);
369 h->SetStats(0);
370 return h;
371}
372
373
374// Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
375
376void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
377 const Double_t er[kNZ][kNR], Double_t &er_value )
378{
379 //
380 // Interpolate table - 2D interpolation
381 //
382 Double_t save_er[10] ;
383
384 Search( kNZ, fgkZList, z, fJLow ) ;
385 Search( kNR, fgkRList, r, fKLow ) ;
386 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
387 if ( fKLow < 0 ) fKLow = 0 ;
388 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
389 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
390
391 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
392 save_er[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
393 }
394 er_value = Interpolate( &fgkZList[fJLow], save_er, order, z ) ;
395
396}
397
398
399Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
400 const Int_t order, const Double_t x )
401{
402 //
403 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
404 //
405
406 Double_t y ;
407 if ( order == 2 ) { // Quadratic Interpolation = 2
408 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
409 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
410 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
411 } else { // Linear Interpolation = 1
412 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
413 }
414
415 return (y);
416
417}
418
419
420void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low )
421{
422 //
423 // Search an ordered table by starting at the most recently used point
424 //
425
426 Long_t middle, high ;
427 Int_t ascend = 0, increment = 1 ;
428
429 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
430
431 if ( low < 0 || low > n-1 ) {
432 low = -1 ; high = n ;
433 } else { // Ordered Search phase
434 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
435 if ( low == n-1 ) return ;
436 high = low + 1 ;
437 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
438 low = high ;
439 increment *= 2 ;
440 high = low + increment ;
441 if ( high > n-1 ) { high = n ; break ; }
442 }
443 } else {
444 if ( low == 0 ) { low = -1 ; return ; }
445 high = low - 1 ;
446 while ( (Int_t)( x < xArray[low] ) == ascend ) {
447 high = low ;
448 increment *= 2 ;
449 if ( increment >= high ) { low = -1 ; break ; }
450 else low = high - increment ;
451 }
452 }
453 }
454
455 while ( (high-low) != 1 ) { // Binary Search Phase
456 middle = ( high + low ) / 2 ;
457 if ( (Int_t)( x >= xArray[middle] ) == ascend )
458 low = middle ;
459 else
460 high = middle ;
461 }
462
463 if ( x == xArray[n-1] ) low = n-2 ;
464 if ( x == xArray[0] ) low = 0 ;
465
466}
467
cf5b0aa0 468
469AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(const AliExternalTrackParam * trackIn, Double_t refX, Int_t dir,TTreeSRedirector *pcstream){
470 //
471 // Fit the track parameters - without and with distortion
472 // 1. Space points in the TPC are simulated along the trajectory
473 // 2. Space points distorted
474 // 3. Fits the non distorted and distroted track to the reference plane at refX
475 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
476 //
477 // trackIn - input track parameters
478 // refX - reference X to fit the track
479 // dir - direction - out=1 or in=-1
480 // pcstream - debug streamer to check the results
481 //
cad404e1 482 // see AliExternalTrackParam.h documentation:
483 // track1.fP[0] - local y (rphi)
484 // track1.fP[1] - z
485 // track1.fP[2] - sinus of local inclination angle
486 // track1.fP[3] - tangent of deep angle
487 // track1.fP[4] - 1/pt
cf5b0aa0 488 AliTPCROC * roc = AliTPCROC::Instance();
489 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
490 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
491 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
492
493 const Double_t kMaxSnp = 0.85;
494 const Double_t kSigmaY=0.1;
495 const Double_t kSigmaZ=0.1;
496 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
497
498 AliExternalTrackParam track(*trackIn); //
499 // generate points
500 AliTrackPointArray pointArray0(npoints0);
501 AliTrackPointArray pointArray1(npoints0);
502 Double_t xyz[3];
503 AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
504 //
505 // simulate the track
506 Int_t npoints=0;
507 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
508 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
509 AliTrackerBase::PropagateTrackTo(&track,radius,kMass,3,kTRUE,kMaxSnp);
510 track.GetXYZ(xyz);
511 AliTrackPoint pIn0; // space point
512 AliTrackPoint pIn1;
513 Int_t roc= (xyz[2]>0)? 0:18;
514 pointArray0.GetPoint(pIn0,npoints);
515 pointArray1.GetPoint(pIn1,npoints);
516 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
517 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
518 DistortPoint(distPoint, roc);
519 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
520 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
521 //
522 track.Rotate(alpha);
523 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
524 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
525 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
526 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
527 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
528 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
529 pointArray0.AddPoint(npoints, &pIn0);
530 pointArray1.AddPoint(npoints, &pIn1);
531 npoints++;
532 if (npoints>=npoints0) break;
533 }
534 //
535 // refit track
536 //
537 AliExternalTrackParam *track0=0;
538 AliExternalTrackParam *track1=0;
539 AliTrackPoint point1,point2,point3;
540 if (dir==1) { //make seed inner
541 pointArray0.GetPoint(point1,1);
542 pointArray0.GetPoint(point2,10);
543 pointArray0.GetPoint(point3,20);
544 }
545 if (dir==-1){ //make seed outer
546 pointArray0.GetPoint(point1,npoints-20);
547 pointArray0.GetPoint(point2,npoints-10);
548 pointArray0.GetPoint(point3,npoints-1);
549 }
550 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
551 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
552
553
554 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
555 Int_t ipoint= (dir>0) ? ipoint: npoints-1-jpoint;
556 //
557 AliTrackPoint pIn0;
558 AliTrackPoint pIn1;
559 pointArray0.GetPoint(pIn0,ipoint);
560 pointArray1.GetPoint(pIn1,ipoint);
561 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
562 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
563 //
564 AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,1,kFALSE,kMaxSnp);
565 AliTrackerBase::PropagateTrackTo(track1,prot1.GetX(),kMass,1,kFALSE,kMaxSnp);
566 track.GetXYZ(xyz);
567 //
568 Double_t pointPos[2]={0,0};
569 Double_t pointCov[3]={0,0,0};
570 pointPos[0]=prot0.GetY();//local y
571 pointPos[1]=prot0.GetZ();//local z
572 pointCov[0]=prot0.GetCov()[3];//simay^2
573 pointCov[1]=prot0.GetCov()[4];//sigmayz
574 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
575 track0->Update(pointPos,pointCov);
576 //
577 pointPos[0]=prot1.GetY();//local y
578 pointPos[1]=prot1.GetZ();//local z
579 pointCov[0]=prot1.GetCov()[3];//simay^2
580 pointCov[1]=prot1.GetCov()[4];//sigmayz
581 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
582 track1->Update(pointPos,pointCov);
583 }
584
585 AliTrackerBase::PropagateTrackTo(track0,refX,kMass,2.,kTRUE,kMaxSnp);
586 track1->Rotate(track0->GetAlpha());
587 AliTrackerBase::PropagateTrackTo(track1,refX,kMass,2.,kFALSE,kMaxSnp);
588
cad404e1 589 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
cf5b0aa0 590 "point0.="<<&pointArray0<< // points
591 "point1.="<<&pointArray1<< // distorted points
592 "trackIn.="<<&track<< // original track
593 "track0.="<<track0<< // fitted track
594 "track1.="<<track1<< // fitted distorted track
595 "\n";
596 delete track0;
597 return track1;
598}
599
600