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