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