]> 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
Adding the correction for the tracking efficiency for the net-protons
[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
b1f0a2a5 68const Double_t AliTPCCorrection::fgkTPCZ0 =249.7; // nominal gating grid position
0116859c 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,
b1f0a2a5 394 const Double_t er[kNZ][kNR], Double_t &erValue ) {
0116859c 395 //
396 // Interpolate table - 2D interpolation
397 //
b1f0a2a5 398 Double_t saveEr[10] ;
0116859c 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++ ) {
b1f0a2a5 408 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
0116859c 409 }
b1f0a2a5 410 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
0116859c 411
412}
413
414
415Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
b1f0a2a5 416 const Int_t order, const Double_t x ) {
0116859c 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
b1f0a2a5 435void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
0116859c 436 //
437 // Search an ordered table by starting at the most recently used point
438 //
439
440 Long_t middle, high ;
441 Int_t ascend = 0, increment = 1 ;
442
443 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
444
445 if ( low < 0 || low > n-1 ) {
446 low = -1 ; high = n ;
447 } else { // Ordered Search phase
448 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
449 if ( low == n-1 ) return ;
450 high = low + 1 ;
451 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
452 low = high ;
453 increment *= 2 ;
454 high = low + increment ;
455 if ( high > n-1 ) { high = n ; break ; }
456 }
457 } else {
458 if ( low == 0 ) { low = -1 ; return ; }
459 high = low - 1 ;
460 while ( (Int_t)( x < xArray[low] ) == ascend ) {
461 high = low ;
462 increment *= 2 ;
463 if ( increment >= high ) { low = -1 ; break ; }
464 else low = high - increment ;
465 }
466 }
467 }
468
469 while ( (high-low) != 1 ) { // Binary Search Phase
470 middle = ( high + low ) / 2 ;
471 if ( (Int_t)( x >= xArray[middle] ) == ascend )
472 low = middle ;
473 else
474 high = middle ;
475 }
476
477 if ( x == xArray[n-1] ) low = n-2 ;
478 if ( x == xArray[0] ) low = 0 ;
479
480}
481
cf5b0aa0 482
b1f0a2a5 483AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
cf5b0aa0 484 //
485 // Fit the track parameters - without and with distortion
486 // 1. Space points in the TPC are simulated along the trajectory
487 // 2. Space points distorted
488 // 3. Fits the non distorted and distroted track to the reference plane at refX
489 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
490 //
491 // trackIn - input track parameters
492 // refX - reference X to fit the track
493 // dir - direction - out=1 or in=-1
494 // pcstream - debug streamer to check the results
495 //
cad404e1 496 // see AliExternalTrackParam.h documentation:
497 // track1.fP[0] - local y (rphi)
498 // track1.fP[1] - z
499 // track1.fP[2] - sinus of local inclination angle
500 // track1.fP[3] - tangent of deep angle
501 // track1.fP[4] - 1/pt
cf5b0aa0 502 AliTPCROC * roc = AliTPCROC::Instance();
503 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
504 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
505 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
506
507 const Double_t kMaxSnp = 0.85;
508 const Double_t kSigmaY=0.1;
509 const Double_t kSigmaZ=0.1;
510 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
511
be67055b 512 AliExternalTrackParam track(trackIn); //
cf5b0aa0 513 // generate points
514 AliTrackPointArray pointArray0(npoints0);
515 AliTrackPointArray pointArray1(npoints0);
516 Double_t xyz[3];
8b63d99c 517 AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
cf5b0aa0 518 //
519 // simulate the track
520 Int_t npoints=0;
521 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
522 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
8b63d99c 523 AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp);
cf5b0aa0 524 track.GetXYZ(xyz);
4486a91f 525 xyz[0]+=gRandom->Gaus(0,0.001);
526 xyz[1]+=gRandom->Gaus(0,0.001);
cf5b0aa0 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);
4486a91f 558 pointArray0.GetPoint(point2,30);
559 pointArray0.GetPoint(point3,60);
cf5b0aa0 560 }
561 if (dir==-1){ //make seed outer
4486a91f 562 pointArray0.GetPoint(point1,npoints-60);
563 pointArray0.GetPoint(point2,npoints-30);
cf5b0aa0 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 //
b322e06a 580 AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
581 AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
8b63d99c 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
b1f0a2a5 688void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
be67055b 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
b322e06a 705 const Double_t kMaxSnp = 0.85;
706 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
707 // const Double_t kB2C=-0.299792458e-3;
be67055b 708 const Int_t kMinEntries=50;
709 Double_t phi,theta, snp, mean,rms, entries;
710 tinput->SetBranchAddress("theta",&theta);
711 tinput->SetBranchAddress("phi", &phi);
712 tinput->SetBranchAddress("snp",&snp);
713 tinput->SetBranchAddress("mean",&mean);
714 tinput->SetBranchAddress("rms",&rms);
715 tinput->SetBranchAddress("entries",&entries);
716 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
717 //
718 Int_t nentries=tinput->GetEntries();
719 Int_t ncorr=corrArray->GetEntries();
720 Double_t corrections[100]; //
721 Double_t tPar[5];
722 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
723 Double_t refX=0;
724 Int_t dir=0;
725 if (dtype==0) {refX=85; dir=-1;}
8b63d99c 726 if (dtype==1) {refX=275; dir=1;}
b322e06a 727 if (dtype==2) {refX=85; dir=-1;}
be67055b 728 //
729 for (Int_t ientry=0; ientry<nentries; ientry+=step){
730 tinput->GetEntry(ientry);
731 tPar[0]=0;
732 tPar[1]=theta*refX;
733 tPar[2]=snp;
734 tPar[3]=theta;
4486a91f 735 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
8b63d99c 736 Double_t bz=AliTrackerBase::GetBz();
4486a91f 737 if (refX>10. && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
738 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
be67055b 739 if (TMath::Abs(snp)>0.251) continue;
740 (*pcstream)<<"fit"<<
8b63d99c 741 "bz="<<bz<< // magnetic filed used
be67055b 742 "dtype="<<dtype<< // detector match type
743 "ptype="<<ptype<< // parameter type
744 "theta="<<theta<< // theta
745 "phi="<<phi<< // phi
746 "snp="<<snp<< // snp
747 "mean="<<mean<< // mean dist value
748 "rms="<<rms<< // rms
749 "entries="<<entries;// number of entries in bin
750 //
751 for (Int_t icorr=0; icorr<ncorr; icorr++) {
752 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
753 corrections[icorr]=0;
754 if (entries>kMinEntries){
755 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
756 AliExternalTrackParam *trackOut = 0;
757 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
758 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
b1f0a2a5 759 if (dtype==0) {refX=85; dir=-1;}
760 if (dtype==1) {refX=275; dir=1;}
761 if (dtype==2) {refX=0; dir=-1;}
762 //
763 AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
764 AliTrackerBase::PropagateTrackToBxByBz(trackOut,refX,kMass,3,kTRUE,kMaxSnp);
b322e06a 765 //
be67055b 766 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
767 delete trackOut;
768 }
769 (*pcstream)<<"fit"<<
770 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
771 }
772 (*pcstream)<<"fit"<<"\n";
773 }
774 delete pcstream;
775}
776
777
778
779
b1f0a2a5 780void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
8b63d99c 781 //
782 // make a distortion map out ou fthe residual histogram
783 // Results are written to the debug streamer - pcstream
784 // Parameters:
785 // his0 - input (4D) residual histogram
786 // pcstream - file to write the tree
787 // run - run number
788 // marian.ivanov@cern.ch
789 const Int_t kMinEntries=50;
790 Int_t nbins1=his0->GetAxis(1)->GetNbins();
791 Int_t first1=his0->GetAxis(1)->GetFirst();
792 Int_t last1 =his0->GetAxis(1)->GetLast();
793 //
794 Double_t bz=AliTrackerBase::GetBz();
795 Int_t idim[4]={0,1,2,3};
796 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
797 //
798 his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
799 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
800 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
801 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
802 Int_t first3 = his1->GetAxis(3)->GetFirst();
803 Int_t last3 = his1->GetAxis(3)->GetLast();
804 //
805
806 for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){ // axis 3 - local angle
807 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
808 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
809 if (ibin3<first3) {
810 his1->GetAxis(3)->SetRangeUser(-1,1);
811 x3=0;
812 }
813 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
814 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
815 Int_t first2 = his3->GetAxis(2)->GetFirst();
816 Int_t last2 = his3->GetAxis(2)->GetLast();
817 //
818 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
819 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
820 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
821 TH1 * hisDelta = his3->Projection(0);
822 //
823 Double_t entries = hisDelta->GetEntries();
824 Double_t mean=0, rms=0;
825 if (entries>kMinEntries){
826 mean = hisDelta->GetMean();
827 rms = hisDelta->GetRMS();
828 }
829 (*pcstream)<<hname<<
830 "run="<<run<<
831 "bz="<<bz<<
832 "theta="<<x1<<
833 "phi="<<x2<<
834 "snp="<<x3<<
835 "entries="<<entries<<
836 "mean="<<mean<<
837 "rms="<<rms<<
838 "\n";
839 delete hisDelta;
840 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
841 }
842 delete his3;
843 }
844 delete his1;
845 }
846}
847
848
849
850
851
ffab0c37 852void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
853 //
854 // Store object in the OCDB
855 // By default the object is stored in the current directory
856 // default comment consit of user name and the date
857 //
858 TString ocdbStorage="";
859 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
860 AliCDBMetaData *metaData= new AliCDBMetaData();
861 metaData->SetObjectClassName("AliTPCCorrection");
862 metaData->SetResponsible("Marian Ivanov");
863 metaData->SetBeamPeriod(1);
864 metaData->SetAliRootVersion("05-25-01"); //root version
865 TString userName=gSystem->GetFromPipe("echo $USER");
866 TString date=gSystem->GetFromPipe("date");
867
868 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
869 if (comment) metaData->SetComment(comment);
870 AliCDBId* id1=NULL;
871 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
872 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
873 gStorage->Put(this, (*id1), metaData);
874}
875