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