1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ////////////////////////////////////////////////////////////////////////////////
18 // AliTPCCorrection class //
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 //
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. //
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 //
33 // The class allows "effective Omega Tau" corrections to be shifted to the //
34 // single distortion classes. //
36 // Note: This class is normally used via the class AliTPCComposedCorrection //
38 // date: 27/04/2010 //
39 // Authors: Magnus Mager, Stefan Rossegger, Jim Thomas //
40 ////////////////////////////////////////////////////////////////////////////////
41 #include "Riostream.h"
46 #include <TTreeStream.h>
49 #include <TTimeStamp.h>
50 #include <AliCDBStorage.h>
52 #include <AliCDBMetaData.h>
57 #include "AliExternalTrackParam.h"
58 #include "AliTrackPointArray.h"
59 #include "TDatabasePDG.h"
60 #include "AliTrackerBase.h"
61 #include "AliTPCROC.h"
62 #include "THnSparse.h"
65 #include "AliTPCTransform.h"
66 #include "AliTPCcalibDB.h"
67 #include "AliTPCExB.h"
68 #include "AliTPCCorrection.h"
69 #include "AliTPCRecoParam.h"
71 #include "AliExternalTrackParam.h"
72 #include "AliTrackPointArray.h"
73 #include "TDatabasePDG.h"
74 #include "AliTrackerBase.h"
75 #include "AliTPCROC.h"
76 #include "THnSparse.h"
78 #include "AliTPCLaserTrack.h"
80 #include "AliTPCCorrection.h"
83 ClassImp(AliTPCCorrection)
85 // FIXME: the following values should come from the database
86 const Double_t AliTPCCorrection::fgkTPCZ0 =249.7; // nominal gating grid position
87 const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
88 const Double_t AliTPCCorrection::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
89 const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
90 const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
91 const Double_t AliTPCCorrection::fgkGG =-70.0; // Gating Grid voltage (volts)
94 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
95 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
96 84.0, 84.5, 85.0, 85.5, 86.0, 87.0, 88.0,
97 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
98 110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
99 130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
100 150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
101 170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
102 190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
103 210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
104 230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 242.0, 244.0, 246.0, 248.0,
105 249.0, 249.5, 250.0, 251.5, 252.0 } ;
107 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = {
108 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
109 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
110 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
111 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
112 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
113 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
114 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
115 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
116 -100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0,
117 -80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0,
118 -60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0,
119 -40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0,
120 -20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0,
121 -1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0,
122 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
123 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0,
124 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0,
125 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0,
126 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0,
127 102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
128 122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
129 142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
130 162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
131 182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
132 202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0,
133 222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0,
134 242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ;
138 AliTPCCorrection::AliTPCCorrection()
139 : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
142 // default constructor
146 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
147 : TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
150 // default constructor, that set the name and title
154 AliTPCCorrection::~AliTPCCorrection() {
156 // virtual destructor
160 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
162 // Corrects the initial coordinates x (cartesian coordinates)
163 // according to the given effect (inherited classes)
164 // roc represents the TPC read out chamber (offline numbering convention)
167 GetCorrection(x,roc,dx);
168 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
171 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
173 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
174 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
175 // roc represents the TPC read out chamber (offline numbering convention)
178 GetCorrection(x,roc,dx);
179 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
182 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
184 // Distorts the initial coordinates x (cartesian coordinates)
185 // according to the given effect (inherited classes)
186 // roc represents the TPC read out chamber (offline numbering convention)
189 GetDistortion(x,roc,dx);
190 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
193 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
195 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
196 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
197 // roc represents the TPC read out chamber (offline numbering convention)
200 GetDistortion(x,roc,dx);
201 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
204 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
206 // This function delivers the correction values dx in respect to the inital coordinates x
207 // roc represents the TPC read out chamber (offline numbering convention)
208 // Note: The dx is overwritten by the inherited effectice class ...
210 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
213 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
215 // This function delivers the distortion values dx in respect to the inital coordinates x
216 // roc represents the TPC read out chamber (offline numbering convention)
218 GetCorrection(x,roc,dx);
219 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
222 void AliTPCCorrection::Init() {
224 // Initialization funtion (not used at the moment)
228 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
234 void AliTPCCorrection::Print(Option_t* /*option*/) const {
236 // Print function to check which correction classes are used
237 // option=="d" prints details regarding the setted magnitude
238 // option=="a" prints the C0 and C1 coefficents for calibration purposes
240 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
243 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
245 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
246 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
251 //SetOmegaTauT1T2(omegaTau, t1, t2);
254 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
256 // Simple plot functionality.
257 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
258 // in respect to position z within the XY plane.
259 // The histogramm has nx times ny entries.
262 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
263 nx,-250.,250.,ny,-250.,250.);
266 Int_t roc=z>0.?0:18; // FIXME
267 for (Int_t iy=1;iy<=ny;++iy) {
268 x[1]=h->GetYaxis()->GetBinCenter(iy);
269 for (Int_t ix=1;ix<=nx;++ix) {
270 x[0]=h->GetXaxis()->GetBinCenter(ix);
271 GetCorrection(x,roc,dx);
272 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
273 if (90.<=r0 && r0<=250.) {
274 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
275 h->SetBinContent(ix,iy,r1-r0);
278 h->SetBinContent(ix,iy,0.);
284 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
286 // Simple plot functionality.
287 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
288 // in respect to position z within the XY plane.
289 // The histogramm has nx times ny entries.
292 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
293 nx,-250.,250.,ny,-250.,250.);
296 Int_t roc=z>0.?0:18; // FIXME
297 for (Int_t iy=1;iy<=ny;++iy) {
298 x[1]=h->GetYaxis()->GetBinCenter(iy);
299 for (Int_t ix=1;ix<=nx;++ix) {
300 x[0]=h->GetXaxis()->GetBinCenter(ix);
301 GetCorrection(x,roc,dx);
302 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
303 if (90.<=r0 && r0<=250.) {
304 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
305 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
307 Float_t dphi=phi1-phi0;
308 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
309 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
311 h->SetBinContent(ix,iy,r0*dphi);
314 h->SetBinContent(ix,iy,0.);
320 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
322 // Simple plot functionality.
323 // Returns a 2d hisogram which represents the corrections in r direction (dr)
324 // in respect to angle phi within the ZR plane.
325 // The histogramm has nx times ny entries.
327 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
328 nz,-250.,250.,nr,85.,250.);
330 for (Int_t ir=1;ir<=nr;++ir) {
331 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
332 x[0]=radius*TMath::Cos(phi);
333 x[1]=radius*TMath::Sin(phi);
334 for (Int_t iz=1;iz<=nz;++iz) {
335 x[2]=h->GetXaxis()->GetBinCenter(iz);
336 Int_t roc=x[2]>0.?0:18; // FIXME
337 GetCorrection(x,roc,dx);
338 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
339 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
340 h->SetBinContent(iz,ir,r1-r0);
347 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
349 // Simple plot functionality.
350 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
351 // in respect to angle phi within the ZR plane.
352 // The histogramm has nx times ny entries.
354 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
355 nz,-250.,250.,nr,85.,250.);
357 for (Int_t iz=1;iz<=nz;++iz) {
358 x[2]=h->GetXaxis()->GetBinCenter(iz);
359 Int_t roc=x[2]>0.?0:18; // FIXME
360 for (Int_t ir=1;ir<=nr;++ir) {
361 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
362 x[0]=radius*TMath::Cos(phi);
363 x[1]=radius*TMath::Sin(phi);
364 GetCorrection(x,roc,dx);
365 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
366 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
367 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
369 Float_t dphi=phi1-phi0;
370 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
371 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
373 h->SetBinContent(iz,ir,r0*dphi);
379 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
380 const char *xlabel,const char *ylabel,const char *zlabel,
381 Int_t nbinsx,Double_t xlow,Double_t xup,
382 Int_t nbinsy,Double_t ylow,Double_t yup) {
384 // Helper function to create a 2d histogramm of given size
390 while (gDirectory->FindObject(hname.Data())) {
397 TH2F *h=new TH2F(hname.Data(),title,
400 h->GetXaxis()->SetTitle(xlabel);
401 h->GetYaxis()->SetTitle(ylabel);
402 h->GetZaxis()->SetTitle(zlabel);
408 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
410 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
411 const Double_t er[kNZ][kNR], Double_t &erValue ) {
413 // Interpolate table - 2D interpolation
415 Double_t saveEr[10] ;
417 Search( kNZ, fgkZList, z, fJLow ) ;
418 Search( kNR, fgkRList, r, fKLow ) ;
419 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
420 if ( fKLow < 0 ) fKLow = 0 ;
421 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
422 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
424 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
425 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
427 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
432 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
433 const Int_t order, const Double_t x ) {
435 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
439 if ( order == 2 ) { // Quadratic Interpolation = 2
440 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
441 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
442 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
443 } else { // Linear Interpolation = 1
444 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
452 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
454 // Search an ordered table by starting at the most recently used point
457 Long_t middle, high ;
458 Int_t ascend = 0, increment = 1 ;
460 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
462 if ( low < 0 || low > n-1 ) {
463 low = -1 ; high = n ;
464 } else { // Ordered Search phase
465 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
466 if ( low == n-1 ) return ;
468 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
471 high = low + increment ;
472 if ( high > n-1 ) { high = n ; break ; }
475 if ( low == 0 ) { low = -1 ; return ; }
477 while ( (Int_t)( x < xArray[low] ) == ascend ) {
480 if ( increment >= high ) { low = -1 ; break ; }
481 else low = high - increment ;
486 while ( (high-low) != 1 ) { // Binary Search Phase
487 middle = ( high + low ) / 2 ;
488 if ( (Int_t)( x >= xArray[middle] ) == ascend )
494 if ( x == xArray[n-1] ) low = n-2 ;
495 if ( x == xArray[0] ) low = 0 ;
499 void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, const TMatrixD &chargeDensity,
500 TMatrixD &arrayErOverEz, const Int_t rows,
501 const Int_t columns, const Int_t iterations ) {
503 // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
505 // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
506 // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
507 // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
508 // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
509 // you wish to solve Laplaces equation however it should not contain random numbers or you will get
510 // random numbers back as a solution.
511 // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
512 // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
513 // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
514 // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
515 // number of points and solves the problem again. This happens several times until the maximum number
516 // of points has been included in the array.
518 // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
519 // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
521 // Original code by Jim Thomas (STAR TPC Collaboration)
524 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
526 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
527 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
528 const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
530 TMatrixD arrayEr(rows,columns) ;
531 TMatrixD arrayEz(rows,columns) ;
533 //Check that number of rows and columns is suitable for a binary expansion
535 if ( !IsPowerOfTwo(rows-1) ) {
536 AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
539 if ( !IsPowerOfTwo(columns-1) ) {
540 AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
544 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
545 // Allow for different size grid spacing in R and Z directions
546 // Use a binary expansion of the size of the matrix to speed up the solution of the problem
548 Int_t iOne = (rows-1)/4 ;
549 Int_t jOne = (columns-1)/4 ;
550 // Solve for N in 2**N, add one.
551 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
553 for ( Int_t count = 0 ; count < loops ; count++ ) {
554 // Loop while the matrix expands & the resolution increases.
556 Float_t tempGridSizeR = gridSizeR * iOne ;
557 Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
558 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
560 // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
561 std::vector<float> coef1(rows) ;
562 std::vector<float> coef2(rows) ;
564 for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
565 Float_t radius = fgkIFCRadius + i*gridSizeR ;
566 coef1[i] = 1.0 + tempGridSizeR/(2*radius);
567 coef2[i] = 1.0 - tempGridSizeR/(2*radius);
570 TMatrixD sumChargeDensity(rows,columns) ;
572 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
573 Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
574 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
575 if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
577 // Add up all enclosed charge density contributions within 1/2 unit in all directions
578 Float_t weight = 0.0 ;
580 sumChargeDensity(i,j) = 0.0 ;
581 for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
582 for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
583 if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
586 // Note that this is cylindrical geometry
587 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
588 sum += weight*radius ;
591 sumChargeDensity(i,j) /= sum ;
593 sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
597 for ( Int_t k = 1 ; k <= iterations; k++ ) {
598 // Solve Poisson's Equation
599 // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
600 // as interations increase.
601 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
602 Float_t overRelaxM1 = overRelax - 1.0 ;
603 Float_t overRelaxtempFourth, overRelaxcoef5 ;
604 overRelaxtempFourth = overRelax * tempFourth ;
605 overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
607 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
608 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
610 arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
611 + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
612 - overRelaxcoef5 * arrayV(i,j)
613 + coef1[i] * arrayV(i+iOne,j)
614 + sumChargeDensity(i,j)
615 ) * overRelaxtempFourth;
619 if ( k == iterations ) {
620 // After full solution is achieved, copy low resolution solution into higher res array
621 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
622 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
625 arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
626 if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
629 arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
630 if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
632 if ( iOne > 1 && jOne > 1 ) {
633 arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
634 if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
635 if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
636 // Note that this leaves a point at the upper left and lower right corners uninitialized.
637 // -> Not a big deal.
646 iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
647 jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
651 // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
652 for ( Int_t j = 0 ; j < columns ; j++ ) {
653 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
654 arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
655 arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
658 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
659 for ( Int_t i = 0 ; i < rows ; i++) {
660 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
661 arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
662 arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
665 for ( Int_t i = 0 ; i < rows ; i++) {
666 // Note: go back and compare to old version of this code. See notes below.
667 // JT Test ... attempt to divide by real Ez not Ez to first order
668 for ( Int_t j = 0 ; j < columns ; j++ ) {
669 arrayEz(i,j) += ezField;
670 // This adds back the overall Z gradient of the field (main E field component)
672 // Warning: (-=) assumes you are using an error potetial without the overall Field included
675 // Integrate Er/Ez from Z to zero
676 for ( Int_t j = 0 ; j < columns ; j++ ) {
677 for ( Int_t i = 0 ; i < rows ; i++ ) {
678 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
679 arrayErOverEz(i,j) = 0.0 ;
680 for ( Int_t k = j ; k < columns ; k++ ) {
681 arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
682 if ( index != 4 ) index = 4; else index = 2 ;
684 if ( index == 4 ) arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
685 if ( index == 2 ) arrayErOverEz(i,j) +=
686 (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
687 -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
688 if ( j == columns-2 ) arrayErOverEz(i,j) =
689 (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
690 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
691 if ( j == columns-1 ) arrayErOverEz(i,j) = 0.0 ;
699 Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
701 // Helperfunction: Check if integer is a power of 2
704 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
705 if ( j == 1 ) return(1) ; // True
710 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
712 // Fit the track parameters - without and with distortion
713 // 1. Space points in the TPC are simulated along the trajectory
714 // 2. Space points distorted
715 // 3. Fits the non distorted and distroted track to the reference plane at refX
716 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
718 // trackIn - input track parameters
719 // refX - reference X to fit the track
720 // dir - direction - out=1 or in=-1
721 // pcstream - debug streamer to check the results
723 // see AliExternalTrackParam.h documentation:
724 // track1.fP[0] - local y (rphi)
726 // track1.fP[2] - sinus of local inclination angle
727 // track1.fP[3] - tangent of deep angle
728 // track1.fP[4] - 1/pt
730 AliTPCROC * roc = AliTPCROC::Instance();
731 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
732 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
733 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
735 const Double_t kMaxSnp = 0.85;
736 const Double_t kSigmaY=0.1;
737 const Double_t kSigmaZ=0.1;
738 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
740 AliExternalTrackParam track(trackIn); //
742 AliTrackPointArray pointArray0(npoints0);
743 AliTrackPointArray pointArray1(npoints0);
745 AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
747 // simulate the track
749 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
750 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
751 AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp);
753 xyz[0]+=gRandom->Gaus(0,0.005);
754 xyz[1]+=gRandom->Gaus(0,0.005);
755 xyz[2]+=gRandom->Gaus(0,0.005);
756 AliTrackPoint pIn0; // space point
758 Int_t sector= (xyz[2]>0)? 0:18;
759 pointArray0.GetPoint(pIn0,npoints);
760 pointArray1.GetPoint(pIn1,npoints);
761 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
762 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
763 DistortPoint(distPoint, sector);
764 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
765 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
768 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
769 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
770 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
771 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
772 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
773 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
774 pointArray0.AddPoint(npoints, &pIn0);
775 pointArray1.AddPoint(npoints, &pIn1);
777 if (npoints>=npoints0) break;
779 if (npoints<npoints0/2) return 0;
783 AliExternalTrackParam *track0=0;
784 AliExternalTrackParam *track1=0;
785 AliTrackPoint point1,point2,point3;
786 if (dir==1) { //make seed inner
787 pointArray0.GetPoint(point1,1);
788 pointArray0.GetPoint(point2,30);
789 pointArray0.GetPoint(point3,60);
791 if (dir==-1){ //make seed outer
792 pointArray0.GetPoint(point1,npoints-60);
793 pointArray0.GetPoint(point2,npoints-30);
794 pointArray0.GetPoint(point3,npoints-1);
796 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
797 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
799 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
800 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
804 pointArray0.GetPoint(pIn0,ipoint);
805 pointArray1.GetPoint(pIn1,ipoint);
806 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
807 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
809 AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
810 AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
811 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
813 Double_t pointPos[2]={0,0};
814 Double_t pointCov[3]={0,0,0};
815 pointPos[0]=prot0.GetY();//local y
816 pointPos[1]=prot0.GetZ();//local z
817 pointCov[0]=prot0.GetCov()[3];//simay^2
818 pointCov[1]=prot0.GetCov()[4];//sigmayz
819 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
820 track0->Update(pointPos,pointCov);
822 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
823 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
824 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
826 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
827 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
828 pointCov[0]=prot1.GetCov()[3];//simay^2
829 pointCov[1]=prot1.GetCov()[4];//sigmayz
830 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
831 track1->Update(pointPos,pointCov);
834 AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
835 track1->Rotate(track0->GetAlpha());
836 track1->PropagateTo(track0->GetX(),AliTrackerBase::GetBz());
838 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
839 "point0.="<<&pointArray0<< // points
840 "point1.="<<&pointArray1<< // distorted points
841 "trackIn.="<<&track<< // original track
842 "track0.="<<track0<< // fitted track
843 "track1.="<<track1<< // fitted distorted track
845 new(&trackIn) AliExternalTrackParam(*track0);
854 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
856 // create the distortion tree on a mesh with granularity given by step
857 // return the tree with distortions at given position
858 // Map is created on the mesh with given step size
860 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
862 for (Double_t x= -250; x<250; x+=step){
863 for (Double_t y= -250; y<250; y+=step){
864 Double_t r = TMath::Sqrt(x*x+y*y);
867 for (Double_t z= -250; z<250; z+=step){
868 Int_t roc=(z>0)?0:18;
872 Double_t phi = TMath::ATan2(y,x);
873 DistortPoint(xyz,roc);
874 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
875 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
876 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
877 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
878 Double_t dx = xyz[0]-x;
879 Double_t dy = xyz[1]-y;
880 Double_t dz = xyz[2]-z;
882 Double_t drphi=(phi1-phi)*r;
883 (*pcstream)<<"distortion"<<
884 "x="<<x<< // original position
889 "x1="<<xyz[0]<< // distorted position
895 "dx="<<dx<< // delta position
905 TFile f(Form("correction%s.root",GetName()));
906 TTree * tree = (TTree*)f.Get("distortion");
907 TTree * tree2= tree->CopyTree("1");
908 tree2->SetName(Form("dist%s",GetName()));
909 tree2->SetDirectory(0);
917 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
920 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
921 // calculates partial distortions
922 // Partial distortion is stored in the resulting tree
923 // Output is storred in the file distortion_<dettype>_<partype>.root
924 // Partial distortion is stored with the name given by correction name
927 // Parameters of function:
928 // input - input tree
929 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex
930 // ppype - parameter type
931 // corrArray - array with partial corrections
932 // step - skipe entries - if 1 all entries processed - it is slow
933 // debug 0 if debug on also space points dumped - it is slow
934 const Double_t kMaxSnp = 0.85;
935 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
936 // const Double_t kB2C=-0.299792458e-3;
937 const Int_t kMinEntries=50;
938 Double_t phi,theta, snp, mean,rms, entries;
939 tinput->SetBranchAddress("theta",&theta);
940 tinput->SetBranchAddress("phi", &phi);
941 tinput->SetBranchAddress("snp",&snp);
942 tinput->SetBranchAddress("mean",&mean);
943 tinput->SetBranchAddress("rms",&rms);
944 tinput->SetBranchAddress("entries",&entries);
945 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
947 Int_t nentries=tinput->GetEntries();
948 Int_t ncorr=corrArray->GetEntries();
949 Double_t corrections[100]={0}; //
951 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
954 if (dtype==0) {refX=85.; dir=-1;}
955 if (dtype==1) {refX=275.; dir=1;}
956 if (dtype==2) {refX=85.; dir=-1;}
957 if (dtype==3) {refX=360.; dir=-1;}
959 for (Int_t ientry=0; ientry<nentries; ientry+=step){
960 tinput->GetEntry(ientry);
961 if (TMath::Abs(snp)>kMaxSnp) continue;
966 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
967 Double_t bz=AliTrackerBase::GetBz();
968 if (refX>10. && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
969 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
970 AliExternalTrackParam track(refX,phi,tPar,cov);
974 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
975 if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature
977 "bz="<<bz<< // magnetic filed used
978 "dtype="<<dtype<< // detector match type
979 "ptype="<<ptype<< // parameter type
980 "theta="<<theta<< // theta
983 "mean="<<mean<< // mean dist value
985 "gx="<<xyz[0]<< // global position at reference
986 "gy="<<xyz[1]<< // global position at reference
987 "gz="<<xyz[2]<< // global position at reference
988 "dRrec="<<dRrec<< // delta Radius in reconstruction
989 "id="<<id<< // track id
990 "entries="<<entries;// number of entries in bin
992 for (Int_t icorr=0; icorr<ncorr; icorr++) {
993 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
994 corrections[icorr]=0;
995 if (entries>kMinEntries){
996 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
997 AliExternalTrackParam *trackOut = 0;
998 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
999 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1000 if (dtype==0) {refX=85.; dir=-1;}
1001 if (dtype==1) {refX=275.; dir=1;}
1002 if (dtype==2) {refX=0; dir=-1;}
1003 if (dtype==3) {refX=360.; dir=-1;}
1006 AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1007 trackOut->Rotate(trackIn.GetAlpha());
1008 trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1010 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1013 corrections[icorr]=0;
1015 if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature
1018 (*pcstream)<<"fit"<<
1019 Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
1020 Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
1021 // for points it is neccessary
1023 (*pcstream)<<"fit"<<"\n";
1030 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1032 // Make a laser fit tree for global minimization
1034 const Double_t cutErrY=0.1;
1035 const Double_t cutErrZ=0.1;
1036 const Double_t kEpsilon=0.00000001;
1041 AliTPCLaserTrack *ltr=0;
1042 AliTPCLaserTrack::LoadTracks();
1043 tree->SetBranchAddress("dY.",&vecdY);
1044 tree->SetBranchAddress("dZ.",&vecdZ);
1045 tree->SetBranchAddress("eY.",&veceY);
1046 tree->SetBranchAddress("eZ.",&veceZ);
1047 tree->SetBranchAddress("LTr.",<r);
1048 Int_t entries= tree->GetEntries();
1049 TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1050 Double_t bz=AliTrackerBase::GetBz();
1053 for (Int_t ientry=0; ientry<entries; ientry++){
1054 tree->GetEntry(ientry);
1055 if (!ltr->GetVecGX()){
1056 ltr->UpdatePoints();
1058 TVectorD * delta= (itype==0)? vecdY:vecdZ;
1059 TVectorD * err= (itype==0)? veceY:veceZ;
1061 for (Int_t irow=0; irow<159; irow++){
1062 Int_t nentries = 1000;
1063 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1064 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1066 Double_t phi =(*ltr->GetVecPhi())[irow];
1067 Double_t theta =ltr->GetTgl();
1068 Double_t mean=delta->GetMatrixArray()[irow];
1069 Double_t gx=0,gy=0,gz=0;
1070 Double_t snp = (*ltr->GetVecP2())[irow];
1071 Double_t rms = 0.1+err->GetMatrixArray()[irow];
1072 gx = (*ltr->GetVecGX())[irow];
1073 gy = (*ltr->GetVecGY())[irow];
1074 gz = (*ltr->GetVecGZ())[irow];
1075 Int_t bundle= ltr->GetBundle();
1078 // get delta R used in reconstruction
1079 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
1080 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1081 const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1082 Double_t xyz0[3]={gx,gy,gz};
1083 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1085 // old ExB correction
1087 if(recoParam&&recoParam->GetUseExBCorrection()) {
1088 Double_t xyz1[3]={gx,gy,gz};
1089 calib->GetExB()->Correct(xyz0,xyz1);
1090 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1093 if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1094 Float_t xyz1[3]={gx,gy,gz};
1095 Int_t sector=(gz>0)?0:18;
1096 correction->CorrectPoint(xyz1, sector);
1097 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1102 (*pcstream)<<"fit"<<
1103 "bz="<<bz<< // magnetic filed used
1104 "dtype="<<dtype<< // detector match type
1105 "ptype="<<itype<< // parameter type
1106 "theta="<<theta<< // theta
1107 "phi="<<phi<< // phi
1108 "snp="<<snp<< // snp
1109 "mean="<<mean<< // mean dist value
1110 "rms="<<rms<< // rms
1111 "gx="<<gx<< // global position
1112 "gy="<<gy<< // global position
1113 "gz="<<gz<< // global position
1114 "dRrec="<<dRrec<< // delta Radius in reconstruction
1115 "id="<<bundle<< //bundle
1116 "entries="<<nentries;// number of entries in bin
1119 Double_t ky = TMath::Tan(TMath::ASin(snp));
1120 Int_t ncorr = corrArray->GetEntries();
1121 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
1122 Double_t phi0 = TMath::ATan2(gy,gx);
1123 Double_t distortions[1000]={0};
1124 Double_t distortionsR[1000]={0};
1125 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1126 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1127 Float_t distPoint[3]={gx,gy,gz};
1128 Int_t sector= (gz>0)? 0:18;
1130 corr->DistortPoint(distPoint, sector);
1132 // Double_t value=distPoint[2]-gz;
1134 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1135 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1136 Double_t drphi= r0*(phi1-phi0);
1137 Double_t dr = r1-r0;
1138 distortions[icorr] = drphi-ky*dr;
1139 distortionsR[icorr] = dr;
1141 (*pcstream)<<"fit"<<
1142 Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
1143 Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
1145 (*pcstream)<<"fit"<<"\n";
1153 void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
1155 // make a distortion map out ou fthe residual histogram
1156 // Results are written to the debug streamer - pcstream
1158 // his0 - input (4D) residual histogram
1159 // pcstream - file to write the tree
1161 // marian.ivanov@cern.ch
1162 const Int_t kMinEntries=50;
1163 Int_t nbins1=his0->GetAxis(1)->GetNbins();
1164 Int_t first1=his0->GetAxis(1)->GetFirst();
1165 Int_t last1 =his0->GetAxis(1)->GetLast();
1167 Double_t bz=AliTrackerBase::GetBz();
1168 Int_t idim[4]={0,1,2,3};
1169 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
1171 his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1172 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1173 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
1174 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
1175 Int_t first3 = his1->GetAxis(3)->GetFirst();
1176 Int_t last3 = his1->GetAxis(3)->GetLast();
1179 for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){ // axis 3 - local angle
1180 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1181 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1183 his1->GetAxis(3)->SetRangeUser(-1,1);
1186 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
1187 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1188 Int_t first2 = his3->GetAxis(2)->GetFirst();
1189 Int_t last2 = his3->GetAxis(2)->GetLast();
1191 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1192 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1193 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1194 TH1 * hisDelta = his3->Projection(0);
1196 Double_t entries = hisDelta->GetEntries();
1197 Double_t mean=0, rms=0;
1198 if (entries>kMinEntries){
1199 mean = hisDelta->GetMean();
1200 rms = hisDelta->GetRMS();
1202 (*pcstream)<<hname<<
1208 "entries="<<entries<<
1213 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1225 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1227 // Store object in the OCDB
1228 // By default the object is stored in the current directory
1229 // default comment consit of user name and the date
1231 TString ocdbStorage="";
1232 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1233 AliCDBMetaData *metaData= new AliCDBMetaData();
1234 metaData->SetObjectClassName("AliTPCCorrection");
1235 metaData->SetResponsible("Marian Ivanov");
1236 metaData->SetBeamPeriod(1);
1237 metaData->SetAliRootVersion("05-25-01"); //root version
1238 TString userName=gSystem->GetFromPipe("echo $USER");
1239 TString date=gSystem->GetFromPipe("date");
1241 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1242 if (comment) metaData->SetComment(comment);
1244 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1245 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1246 gStorage->Put(this, (*id1), metaData);