Fix coding rule violations
[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>
45
46#include "AliTPCCorrection.h"
47
48// FIXME: the following values should come from the database
49const Double_t AliTPCCorrection::fgkTPC_Z0 =249.7; // nominal gating grid position
50const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
51const Double_t AliTPCCorrection::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
52const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
53const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
54const Double_t AliTPCCorrection::fgkGG =-70.0; // Gating Grid voltage (volts)
55
56
57// FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
58const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
5984.0, 84.5, 85.0, 85.5, 86.0, 87.0, 88.0,
6090.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
61110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
62130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
63150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
64170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
65190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
66210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
67230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 242.0, 244.0, 246.0, 248.0,
68249.0, 249.5, 250.0, 251.5, 252.0 } ;
69
70const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = {
71-249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
72-240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
73-220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
74-200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
75-180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
76-160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
77-140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
78-120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
79-100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0,
80-80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0,
81-60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0,
82-40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0,
83-20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0,
84-1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0,
85 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
86 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0,
87 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0,
88 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0,
89 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0,
90102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
91122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
92142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
93162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
94182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
95202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0,
96222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0,
97242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ;
98
99
100
101AliTPCCorrection::AliTPCCorrection()
102 : TNamed("correction_unity","unity"),fJLow(0),fKLow(0)
103{
104 //
105 // default constructor
106 //
107}
108
109AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
110 : TNamed(name,title),fJLow(0),fKLow(0)
111{
112 //
113 // default constructor, that set the name and title
114 //
115}
116
117AliTPCCorrection::~AliTPCCorrection() {
118 //
119 // virtual destructor
120 //
121}
122
123void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
124 //
125 // Corrects the initial coordinates x (cartesian coordinates)
126 // according to the given effect (inherited classes)
127 // roc represents the TPC read out chamber (offline numbering convention)
128 //
129 Float_t dx[3];
130 GetCorrection(x,roc,dx);
131 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
132}
133
134void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
135 //
136 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
137 // (distorted) coordinates in xp. The distortion is set 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) xp[j]=x[j]+dx[j];
143}
144
145void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
146 //
147 // Distorts the initial coordinates x (cartesian coordinates)
148 // 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 GetDistortion(x,roc,dx);
153 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
154}
155
156void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
157 //
158 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
159 // (distorted) coordinates in xp. The distortion is set 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) xp[j]=x[j]+dx[j];
165}
166
167void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
168 //
169 // This function delivers the correction values dx in respect to the inital coordinates x
170 // roc represents the TPC read out chamber (offline numbering convention)
171 // Note: The dx is overwritten by the inherited effectice class ...
172 //
173 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
174}
175
176void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
177 //
178 // This function delivers the distortion values dx in respect to the inital coordinates x
179 // roc represents the TPC read out chamber (offline numbering convention)
180 //
181 GetCorrection(x,roc,dx);
182 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
183}
184
185void AliTPCCorrection::Init() {
186 //
187 // Initialization funtion (not used at the moment)
188 //
189}
190
191void AliTPCCorrection::Print(Option_t* /*option*/) const {
192 //
193 // Print function to check which correction classes are used
194 // option=="d" prints details regarding the setted magnitude
195 // option=="a" prints the C0 and C1 coefficents for calibration purposes
196 //
197 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
198}
199
200void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t /*t1*/,Float_t /*t2*/) {
201 //
202 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
203 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
204 // calibration run
205 //
206 // SetOmegaTauT1T2(omegaTau, t1, t2);
207}
208
209TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
210 //
211 // Simple plot functionality.
212 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
213 // in respect to position z within the XY plane.
214 // The histogramm has nx times ny entries.
215 //
216
217 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
218 nx,-250.,250.,ny,-250.,250.);
219 Float_t x[3],dx[3];
220 x[2]=z;
221 Int_t roc=z>0.?0:18; // FIXME
222 for (Int_t iy=1;iy<=ny;++iy) {
223 x[1]=h->GetYaxis()->GetBinCenter(iy);
224 for (Int_t ix=1;ix<=nx;++ix) {
225 x[0]=h->GetXaxis()->GetBinCenter(ix);
226 GetCorrection(x,roc,dx);
227 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
228 if (90.<=r0 && r0<=250.) {
229 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
230 h->SetBinContent(ix,iy,r1-r0);
231 }
232 else
233 h->SetBinContent(ix,iy,0.);
234 }
235 }
236 return h;
237}
238
239TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
240 //
241 // Simple plot functionality.
242 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
243 // in respect to position z within the XY plane.
244 // The histogramm has nx times ny entries.
245 //
246
247 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
248 nx,-250.,250.,ny,-250.,250.);
249 Float_t x[3],dx[3];
250 x[2]=z;
251 Int_t roc=z>0.?0:18; // FIXME
252 for (Int_t iy=1;iy<=ny;++iy) {
253 x[1]=h->GetYaxis()->GetBinCenter(iy);
254 for (Int_t ix=1;ix<=nx;++ix) {
255 x[0]=h->GetXaxis()->GetBinCenter(ix);
256 GetCorrection(x,roc,dx);
257 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
258 if (90.<=r0 && r0<=250.) {
259 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
260 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
261
262 Float_t dphi=phi1-phi0;
263 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
264 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
265
266 h->SetBinContent(ix,iy,r0*dphi);
267 }
268 else
269 h->SetBinContent(ix,iy,0.);
270 }
271 }
272 return h;
273}
274
275TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
276 //
277 // Simple plot functionality.
278 // Returns a 2d hisogram which represents the corrections in r direction (dr)
279 // in respect to angle phi within the ZR plane.
280 // The histogramm has nx times ny entries.
281 //
282 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
283 nz,-250.,250.,nr,85.,250.);
284 Float_t x[3],dx[3];
285 for (Int_t ir=1;ir<=nr;++ir) {
286 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
287 x[0]=radius*TMath::Cos(phi);
288 x[1]=radius*TMath::Sin(phi);
289 for (Int_t iz=1;iz<=nz;++iz) {
290 x[2]=h->GetXaxis()->GetBinCenter(iz);
291 Int_t roc=x[2]>0.?0:18; // FIXME
292 GetCorrection(x,roc,dx);
293 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
294 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
295 h->SetBinContent(iz,ir,r1-r0);
296 }
297 }
298 printf("SDF\n");
299 return h;
300
301}
302
303TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
304 //
305 // Simple plot functionality.
306 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
307 // in respect to angle phi within the ZR plane.
308 // The histogramm has nx times ny entries.
309 //
310 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
311 nz,-250.,250.,nr,85.,250.);
312 Float_t x[3],dx[3];
313 for (Int_t iz=1;iz<=nz;++iz) {
314 x[2]=h->GetXaxis()->GetBinCenter(iz);
315 Int_t roc=x[2]>0.?0:18; // FIXME
316 for (Int_t ir=1;ir<=nr;++ir) {
317 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
318 x[0]=radius*TMath::Cos(phi);
319 x[1]=radius*TMath::Sin(phi);
320 GetCorrection(x,roc,dx);
321 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
322 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
323 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
324
325 Float_t dphi=phi1-phi0;
326 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
327 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
328
329 h->SetBinContent(iz,ir,r0*dphi);
330 }
331 }
332 return h;
333}
334
335TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
336 const char *xlabel,const char *ylabel,const char *zlabel,
337 Int_t nbinsx,Double_t xlow,Double_t xup,
338 Int_t nbinsy,Double_t ylow,Double_t yup) {
339 //
340 // Helper function to create a 2d histogramm of given size
341 //
342
343 TString hname=name;
344 Int_t i=0;
345 if (gDirectory) {
346 while (gDirectory->FindObject(hname.Data())) {
347 hname =name;
348 hname+="_";
349 hname+=i;
350 ++i;
351 }
352 }
353 TH2F *h=new TH2F(hname.Data(),title,
354 nbinsx,xlow,xup,
355 nbinsy,ylow,yup);
356 h->GetXaxis()->SetTitle(xlabel);
357 h->GetYaxis()->SetTitle(ylabel);
358 h->GetZaxis()->SetTitle(zlabel);
359 h->SetStats(0);
360 return h;
361}
362
363
364// Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
365
366void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
367 const Double_t er[kNZ][kNR], Double_t &er_value )
368{
369 //
370 // Interpolate table - 2D interpolation
371 //
372 Double_t save_er[10] ;
373
374 Search( kNZ, fgkZList, z, fJLow ) ;
375 Search( kNR, fgkRList, r, fKLow ) ;
376 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
377 if ( fKLow < 0 ) fKLow = 0 ;
378 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
379 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
380
381 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
382 save_er[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
383 }
384 er_value = Interpolate( &fgkZList[fJLow], save_er, order, z ) ;
385
386}
387
388
389Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
390 const Int_t order, const Double_t x )
391{
392 //
393 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
394 //
395
396 Double_t y ;
397 if ( order == 2 ) { // Quadratic Interpolation = 2
398 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
399 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
400 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
401 } else { // Linear Interpolation = 1
402 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
403 }
404
405 return (y);
406
407}
408
409
410void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low )
411{
412 //
413 // Search an ordered table by starting at the most recently used point
414 //
415
416 Long_t middle, high ;
417 Int_t ascend = 0, increment = 1 ;
418
419 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
420
421 if ( low < 0 || low > n-1 ) {
422 low = -1 ; high = n ;
423 } else { // Ordered Search phase
424 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
425 if ( low == n-1 ) return ;
426 high = low + 1 ;
427 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
428 low = high ;
429 increment *= 2 ;
430 high = low + increment ;
431 if ( high > n-1 ) { high = n ; break ; }
432 }
433 } else {
434 if ( low == 0 ) { low = -1 ; return ; }
435 high = low - 1 ;
436 while ( (Int_t)( x < xArray[low] ) == ascend ) {
437 high = low ;
438 increment *= 2 ;
439 if ( increment >= high ) { low = -1 ; break ; }
440 else low = high - increment ;
441 }
442 }
443 }
444
445 while ( (high-low) != 1 ) { // Binary Search Phase
446 middle = ( high + low ) / 2 ;
447 if ( (Int_t)( x >= xArray[middle] ) == ascend )
448 low = middle ;
449 else
450 high = middle ;
451 }
452
453 if ( x == xArray[n-1] ) low = n-2 ;
454 if ( x == xArray[0] ) low = 0 ;
455
456}
457
458ClassImp(AliTPCCorrection)