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