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