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"
79 #include "AliESDVertex.h"
80 #include "AliVertexerTracks.h"
81 #include "TDatabasePDG.h"
84 #include "AliTPCCorrection.h"
87 ClassImp(AliTPCCorrection)
90 TObjArray *AliTPCCorrection::fgVisualCorrection=0;
91 // instance of correction for visualization
94 // FIXME: the following values should come from the database
95 const Double_t AliTPCCorrection::fgkTPCZ0 =249.7; // nominal gating grid position
96 const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
97 const Double_t AliTPCCorrection::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
98 const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
99 const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
100 const Double_t AliTPCCorrection::fgkGG =-70.0; // Gating Grid voltage (volts)
104 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
105 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
106 84.0, 84.5, 85.0, 85.5, 86.0, 87.0, 88.0,
107 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
108 110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
109 130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
110 150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
111 170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
112 190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
113 210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
114 230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 242.0, 244.0, 246.0, 248.0,
115 249.0, 249.5, 250.0, 251.5, 252.0 } ;
117 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = {
118 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
119 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
120 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
121 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
122 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
123 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
124 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
125 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
126 -100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0,
127 -80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0,
128 -60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0,
129 -40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0,
130 -20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0,
131 -1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0,
132 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
133 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0,
134 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0,
135 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0,
136 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0,
137 102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
138 122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
139 142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
140 162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
141 182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
142 202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0,
143 222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0,
144 242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ;
148 AliTPCCorrection::AliTPCCorrection()
149 : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
152 // default constructor
154 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
157 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
158 : TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
161 // default constructor, that set the name and title
163 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
166 AliTPCCorrection::~AliTPCCorrection() {
168 // virtual destructor
172 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
174 // Corrects the initial coordinates x (cartesian coordinates)
175 // according to the given effect (inherited classes)
176 // roc represents the TPC read out chamber (offline numbering convention)
179 GetCorrection(x,roc,dx);
180 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
183 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
185 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
186 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
187 // roc represents the TPC read out chamber (offline numbering convention)
190 GetCorrection(x,roc,dx);
191 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
194 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
196 // Distorts the initial coordinates x (cartesian coordinates)
197 // according to the given effect (inherited classes)
198 // roc represents the TPC read out chamber (offline numbering convention)
201 GetDistortion(x,roc,dx);
202 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
205 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
207 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
208 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
209 // roc represents the TPC read out chamber (offline numbering convention)
212 GetDistortion(x,roc,dx);
213 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
216 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
218 // This function delivers the correction values dx in respect to the inital coordinates x
219 // roc represents the TPC read out chamber (offline numbering convention)
220 // Note: The dx is overwritten by the inherited effectice class ...
222 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
225 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
227 // This function delivers the distortion values dx in respect to the inital coordinates x
228 // roc represents the TPC read out chamber (offline numbering convention)
230 GetCorrection(x,roc,dx);
231 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
234 void AliTPCCorrection::Init() {
236 // Initialization funtion (not used at the moment)
240 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
246 void AliTPCCorrection::Print(Option_t* /*option*/) const {
248 // Print function to check which correction classes are used
249 // option=="d" prints details regarding the setted magnitude
250 // option=="a" prints the C0 and C1 coefficents for calibration purposes
252 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
255 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
257 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
258 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
263 //SetOmegaTauT1T2(omegaTau, t1, t2);
266 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
268 // Simple plot functionality.
269 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
270 // in respect to position z within the XY plane.
271 // The histogramm has nx times ny entries.
274 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
275 nx,-250.,250.,ny,-250.,250.);
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 r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
287 h->SetBinContent(ix,iy,r1-r0);
290 h->SetBinContent(ix,iy,0.);
296 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
298 // Simple plot functionality.
299 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
300 // in respect to position z within the XY plane.
301 // The histogramm has nx times ny entries.
304 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
305 nx,-250.,250.,ny,-250.,250.);
308 Int_t roc=z>0.?0:18; // FIXME
309 for (Int_t iy=1;iy<=ny;++iy) {
310 x[1]=h->GetYaxis()->GetBinCenter(iy);
311 for (Int_t ix=1;ix<=nx;++ix) {
312 x[0]=h->GetXaxis()->GetBinCenter(ix);
313 GetCorrection(x,roc,dx);
314 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
315 if (90.<=r0 && r0<=250.) {
316 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
317 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
319 Float_t dphi=phi1-phi0;
320 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
321 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
323 h->SetBinContent(ix,iy,r0*dphi);
326 h->SetBinContent(ix,iy,0.);
332 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
334 // Simple plot functionality.
335 // Returns a 2d hisogram which represents the corrections in r direction (dr)
336 // in respect to angle phi within the ZR plane.
337 // The histogramm has nx times ny entries.
339 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
340 nz,-250.,250.,nr,85.,250.);
342 for (Int_t ir=1;ir<=nr;++ir) {
343 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
344 x[0]=radius*TMath::Cos(phi);
345 x[1]=radius*TMath::Sin(phi);
346 for (Int_t iz=1;iz<=nz;++iz) {
347 x[2]=h->GetXaxis()->GetBinCenter(iz);
348 Int_t roc=x[2]>0.?0:18; // FIXME
349 GetCorrection(x,roc,dx);
350 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
351 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
352 h->SetBinContent(iz,ir,r1-r0);
359 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
361 // Simple plot functionality.
362 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
363 // in respect to angle phi within the ZR plane.
364 // The histogramm has nx times ny entries.
366 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
367 nz,-250.,250.,nr,85.,250.);
369 for (Int_t iz=1;iz<=nz;++iz) {
370 x[2]=h->GetXaxis()->GetBinCenter(iz);
371 Int_t roc=x[2]>0.?0:18; // FIXME
372 for (Int_t ir=1;ir<=nr;++ir) {
373 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
374 x[0]=radius*TMath::Cos(phi);
375 x[1]=radius*TMath::Sin(phi);
376 GetCorrection(x,roc,dx);
377 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
378 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
379 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
381 Float_t dphi=phi1-phi0;
382 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
383 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
385 h->SetBinContent(iz,ir,r0*dphi);
391 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
392 const char *xlabel,const char *ylabel,const char *zlabel,
393 Int_t nbinsx,Double_t xlow,Double_t xup,
394 Int_t nbinsy,Double_t ylow,Double_t yup) {
396 // Helper function to create a 2d histogramm of given size
402 while (gDirectory->FindObject(hname.Data())) {
409 TH2F *h=new TH2F(hname.Data(),title,
412 h->GetXaxis()->SetTitle(xlabel);
413 h->GetYaxis()->SetTitle(ylabel);
414 h->GetZaxis()->SetTitle(zlabel);
420 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
422 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
423 const Double_t er[kNZ][kNR], Double_t &erValue ) {
425 // Interpolate table - 2D interpolation
427 Double_t saveEr[10] ;
429 Search( kNZ, fgkZList, z, fJLow ) ;
430 Search( kNR, fgkRList, r, fKLow ) ;
431 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
432 if ( fKLow < 0 ) fKLow = 0 ;
433 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
434 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
436 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
437 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
439 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
444 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
445 const Int_t order, const Double_t x ) {
447 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
451 if ( order == 2 ) { // Quadratic Interpolation = 2
452 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
453 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
454 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
455 } else { // Linear Interpolation = 1
456 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
464 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
466 // Search an ordered table by starting at the most recently used point
469 Long_t middle, high ;
470 Int_t ascend = 0, increment = 1 ;
472 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
474 if ( low < 0 || low > n-1 ) {
475 low = -1 ; high = n ;
476 } else { // Ordered Search phase
477 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
478 if ( low == n-1 ) return ;
480 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
483 high = low + increment ;
484 if ( high > n-1 ) { high = n ; break ; }
487 if ( low == 0 ) { low = -1 ; return ; }
489 while ( (Int_t)( x < xArray[low] ) == ascend ) {
492 if ( increment >= high ) { low = -1 ; break ; }
493 else low = high - increment ;
498 while ( (high-low) != 1 ) { // Binary Search Phase
499 middle = ( high + low ) / 2 ;
500 if ( (Int_t)( x >= xArray[middle] ) == ascend )
506 if ( x == xArray[n-1] ) low = n-2 ;
507 if ( x == xArray[0] ) low = 0 ;
511 void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, const TMatrixD &chargeDensity,
512 TMatrixD &arrayErOverEz, const Int_t rows,
513 const Int_t columns, const Int_t iterations ) {
515 // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
517 // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
518 // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
519 // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
520 // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
521 // you wish to solve Laplaces equation however it should not contain random numbers or you will get
522 // random numbers back as a solution.
523 // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
524 // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
525 // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
526 // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
527 // number of points and solves the problem again. This happens several times until the maximum number
528 // of points has been included in the array.
530 // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
531 // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
533 // Original code by Jim Thomas (STAR TPC Collaboration)
536 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
538 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
539 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
540 const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
542 TMatrixD arrayEr(rows,columns) ;
543 TMatrixD arrayEz(rows,columns) ;
545 //Check that number of rows and columns is suitable for a binary expansion
547 if ( !IsPowerOfTwo(rows-1) ) {
548 AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
551 if ( !IsPowerOfTwo(columns-1) ) {
552 AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
556 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
557 // Allow for different size grid spacing in R and Z directions
558 // Use a binary expansion of the size of the matrix to speed up the solution of the problem
560 Int_t iOne = (rows-1)/4 ;
561 Int_t jOne = (columns-1)/4 ;
562 // Solve for N in 2**N, add one.
563 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
565 for ( Int_t count = 0 ; count < loops ; count++ ) {
566 // Loop while the matrix expands & the resolution increases.
568 Float_t tempGridSizeR = gridSizeR * iOne ;
569 Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
570 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
572 // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
573 std::vector<float> coef1(rows) ;
574 std::vector<float> coef2(rows) ;
576 for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
577 Float_t radius = fgkIFCRadius + i*gridSizeR ;
578 coef1[i] = 1.0 + tempGridSizeR/(2*radius);
579 coef2[i] = 1.0 - tempGridSizeR/(2*radius);
582 TMatrixD sumChargeDensity(rows,columns) ;
584 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
585 Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
586 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
587 if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
589 // Add up all enclosed charge density contributions within 1/2 unit in all directions
590 Float_t weight = 0.0 ;
592 sumChargeDensity(i,j) = 0.0 ;
593 for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
594 for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
595 if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
598 // Note that this is cylindrical geometry
599 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
600 sum += weight*radius ;
603 sumChargeDensity(i,j) /= sum ;
605 sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
609 for ( Int_t k = 1 ; k <= iterations; k++ ) {
610 // Solve Poisson's Equation
611 // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
612 // as interations increase.
613 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
614 Float_t overRelaxM1 = overRelax - 1.0 ;
615 Float_t overRelaxtempFourth, overRelaxcoef5 ;
616 overRelaxtempFourth = overRelax * tempFourth ;
617 overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
619 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
620 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
622 arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
623 + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
624 - overRelaxcoef5 * arrayV(i,j)
625 + coef1[i] * arrayV(i+iOne,j)
626 + sumChargeDensity(i,j)
627 ) * overRelaxtempFourth;
631 if ( k == iterations ) {
632 // After full solution is achieved, copy low resolution solution into higher res array
633 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
634 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
637 arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
638 if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
641 arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
642 if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
644 if ( iOne > 1 && jOne > 1 ) {
645 arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
646 if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
647 if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
648 // Note that this leaves a point at the upper left and lower right corners uninitialized.
649 // -> Not a big deal.
658 iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
659 jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
663 // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
664 for ( Int_t j = 0 ; j < columns ; j++ ) {
665 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
666 arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
667 arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
670 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
671 for ( Int_t i = 0 ; i < rows ; i++) {
672 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
673 arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
674 arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
677 for ( Int_t i = 0 ; i < rows ; i++) {
678 // Note: go back and compare to old version of this code. See notes below.
679 // JT Test ... attempt to divide by real Ez not Ez to first order
680 for ( Int_t j = 0 ; j < columns ; j++ ) {
681 arrayEz(i,j) += ezField;
682 // This adds back the overall Z gradient of the field (main E field component)
684 // Warning: (-=) assumes you are using an error potetial without the overall Field included
687 // Integrate Er/Ez from Z to zero
688 for ( Int_t j = 0 ; j < columns ; j++ ) {
689 for ( Int_t i = 0 ; i < rows ; i++ ) {
690 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
691 arrayErOverEz(i,j) = 0.0 ;
692 for ( Int_t k = j ; k < columns ; k++ ) {
693 arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
694 if ( index != 4 ) index = 4; else index = 2 ;
696 if ( index == 4 ) arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
697 if ( index == 2 ) arrayErOverEz(i,j) +=
698 (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
699 -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
700 if ( j == columns-2 ) arrayErOverEz(i,j) =
701 (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
702 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
703 if ( j == columns-1 ) arrayErOverEz(i,j) = 0.0 ;
711 Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
713 // Helperfunction: Check if integer is a power of 2
716 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
717 if ( j == 1 ) return(1) ; // True
722 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
724 // Fit the track parameters - without and with distortion
725 // 1. Space points in the TPC are simulated along the trajectory
726 // 2. Space points distorted
727 // 3. Fits the non distorted and distroted track to the reference plane at refX
728 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
730 // trackIn - input track parameters
731 // refX - reference X to fit the track
732 // dir - direction - out=1 or in=-1
733 // pcstream - debug streamer to check the results
735 // see AliExternalTrackParam.h documentation:
736 // track1.fP[0] - local y (rphi)
738 // track1.fP[2] - sinus of local inclination angle
739 // track1.fP[3] - tangent of deep angle
740 // track1.fP[4] - 1/pt
742 AliTPCROC * roc = AliTPCROC::Instance();
743 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
744 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
745 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
746 const Double_t kMaxSnp = 0.85;
747 const Double_t kSigmaY=0.1;
748 const Double_t kSigmaZ=0.1;
749 const Double_t kMaxR=500;
750 const Double_t kMaxZ=500;
751 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
755 AliExternalTrackParam track(trackIn); //
757 AliTrackPointArray pointArray0(npoints0);
758 AliTrackPointArray pointArray1(npoints0);
760 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
762 // simulate the track
764 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
765 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
766 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
768 xyz[0]+=gRandom->Gaus(0,0.00005);
769 xyz[1]+=gRandom->Gaus(0,0.00005);
770 xyz[2]+=gRandom->Gaus(0,0.00005);
771 if (TMath::Abs(track.GetZ())>kMaxZ) break;
772 if (TMath::Abs(track.GetX())>kMaxR) break;
773 AliTrackPoint pIn0; // space point
775 Int_t sector= (xyz[2]>0)? 0:18;
776 pointArray0.GetPoint(pIn0,npoints);
777 pointArray1.GetPoint(pIn1,npoints);
778 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
779 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
780 DistortPoint(distPoint, sector);
781 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
782 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
785 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
786 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
787 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
788 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
789 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
790 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
791 pointArray0.AddPoint(npoints, &pIn0);
792 pointArray1.AddPoint(npoints, &pIn1);
794 if (npoints>=npoints0) break;
796 if (npoints<npoints0/2) return 0;
800 AliExternalTrackParam *track0=0;
801 AliExternalTrackParam *track1=0;
802 AliTrackPoint point1,point2,point3;
803 if (dir==1) { //make seed inner
804 pointArray0.GetPoint(point1,1);
805 pointArray0.GetPoint(point2,30);
806 pointArray0.GetPoint(point3,60);
808 if (dir==-1){ //make seed outer
809 pointArray0.GetPoint(point1,npoints-60);
810 pointArray0.GetPoint(point2,npoints-30);
811 pointArray0.GetPoint(point3,npoints-1);
813 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
814 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
816 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
817 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
821 pointArray0.GetPoint(pIn0,ipoint);
822 pointArray1.GetPoint(pIn1,ipoint);
823 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
824 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
826 if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
827 if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
828 if (TMath::Abs(track0->GetZ())>kMaxZ) break;
829 if (TMath::Abs(track0->GetX())>kMaxR) break;
830 if (TMath::Abs(track1->GetZ())>kMaxZ) break;
831 if (TMath::Abs(track1->GetX())>kMaxR) break;
833 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
835 Double_t pointPos[2]={0,0};
836 Double_t pointCov[3]={0,0,0};
837 pointPos[0]=prot0.GetY();//local y
838 pointPos[1]=prot0.GetZ();//local z
839 pointCov[0]=prot0.GetCov()[3];//simay^2
840 pointCov[1]=prot0.GetCov()[4];//sigmayz
841 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
842 if (!track0->Update(pointPos,pointCov)) break;
844 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
845 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
846 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
848 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
849 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
850 pointCov[0]=prot1.GetCov()[3];//simay^2
851 pointCov[1]=prot1.GetCov()[4];//sigmayz
852 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
853 if (!track1->Update(pointPos,pointCov)) break;
857 if (npoints2<npoints) return 0;
858 AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
859 track1->Rotate(track0->GetAlpha());
860 AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
862 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
863 "point0.="<<&pointArray0<< // points
864 "point1.="<<&pointArray1<< // distorted points
865 "trackIn.="<<&track<< // original track
866 "track0.="<<track0<< // fitted track
867 "track1.="<<track1<< // fitted distorted track
869 new(&trackIn) AliExternalTrackParam(*track0);
878 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
880 // create the distortion tree on a mesh with granularity given by step
881 // return the tree with distortions at given position
882 // Map is created on the mesh with given step size
884 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
886 for (Double_t x= -250; x<250; x+=step){
887 for (Double_t y= -250; y<250; y+=step){
888 Double_t r = TMath::Sqrt(x*x+y*y);
891 for (Double_t z= -250; z<250; z+=step){
892 Int_t roc=(z>0)?0:18;
896 Double_t phi = TMath::ATan2(y,x);
897 DistortPoint(xyz,roc);
898 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
899 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
900 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
901 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
902 Double_t dx = xyz[0]-x;
903 Double_t dy = xyz[1]-y;
904 Double_t dz = xyz[2]-z;
906 Double_t drphi=(phi1-phi)*r;
907 (*pcstream)<<"distortion"<<
908 "x="<<x<< // original position
913 "x1="<<xyz[0]<< // distorted position
919 "dx="<<dx<< // delta position
929 TFile f(Form("correction%s.root",GetName()));
930 TTree * tree = (TTree*)f.Get("distortion");
931 TTree * tree2= tree->CopyTree("1");
932 tree2->SetName(Form("dist%s",GetName()));
933 tree2->SetDirectory(0);
941 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
944 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
945 // calculates partial distortions
946 // Partial distortion is stored in the resulting tree
947 // Output is storred in the file distortion_<dettype>_<partype>.root
948 // Partial distortion is stored with the name given by correction name
951 // Parameters of function:
952 // input - input tree
953 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex
954 // ppype - parameter type
955 // corrArray - array with partial corrections
956 // step - skipe entries - if 1 all entries processed - it is slow
957 // debug 0 if debug on also space points dumped - it is slow
958 const Double_t kMaxSnp = 0.85;
959 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
960 // const Double_t kB2C=-0.299792458e-3;
961 const Int_t kMinEntries=50;
962 Double_t phi,theta, snp, mean,rms, entries;
963 tinput->SetBranchAddress("theta",&theta);
964 tinput->SetBranchAddress("phi", &phi);
965 tinput->SetBranchAddress("snp",&snp);
966 tinput->SetBranchAddress("mean",&mean);
967 tinput->SetBranchAddress("rms",&rms);
968 tinput->SetBranchAddress("entries",&entries);
969 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
971 Int_t nentries=tinput->GetEntries();
972 Int_t ncorr=corrArray->GetEntries();
973 Double_t corrections[100]={0}; //
975 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
978 if (dtype==0) {refX=85.; dir=-1;}
979 if (dtype==1) {refX=275.; dir=1;}
980 if (dtype==2) {refX=85.; dir=-1;}
981 if (dtype==3) {refX=360.; dir=-1;}
983 for (Int_t ientry=0; ientry<nentries; ientry+=step){
984 tinput->GetEntry(ientry);
985 if (TMath::Abs(snp)>kMaxSnp) continue;
990 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
991 Double_t bz=AliTrackerBase::GetBz();
992 if (refX>10. && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
993 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
994 AliExternalTrackParam track(refX,phi,tPar,cov);
998 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
999 if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature
1000 (*pcstream)<<"fit"<<
1001 "bz="<<bz<< // magnetic filed used
1002 "dtype="<<dtype<< // detector match type
1003 "ptype="<<ptype<< // parameter type
1004 "theta="<<theta<< // theta
1005 "phi="<<phi<< // phi
1006 "snp="<<snp<< // snp
1007 "mean="<<mean<< // mean dist value
1008 "rms="<<rms<< // rms
1009 "gx="<<xyz[0]<< // global position at reference
1010 "gy="<<xyz[1]<< // global position at reference
1011 "gz="<<xyz[2]<< // global position at reference
1012 "dRrec="<<dRrec<< // delta Radius in reconstruction
1013 "id="<<id<< // track id
1014 "entries="<<entries;// number of entries in bin
1016 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1017 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1018 corrections[icorr]=0;
1019 if (entries>kMinEntries){
1020 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1021 AliExternalTrackParam *trackOut = 0;
1022 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1023 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1024 if (dtype==0) {refX=85.; dir=-1;}
1025 if (dtype==1) {refX=275.; dir=1;}
1026 if (dtype==2) {refX=0; dir=-1;}
1027 if (dtype==3) {refX=360.; dir=-1;}
1030 AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1031 trackOut->Rotate(trackIn.GetAlpha());
1032 trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1034 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1037 corrections[icorr]=0;
1039 if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature
1042 (*pcstream)<<"fit"<<
1043 Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
1044 Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
1045 // for points it is neccessary
1047 (*pcstream)<<"fit"<<"\n";
1054 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1056 // Make a laser fit tree for global minimization
1058 const Double_t cutErrY=0.1;
1059 const Double_t cutErrZ=0.1;
1060 const Double_t kEpsilon=0.00000001;
1065 AliTPCLaserTrack *ltr=0;
1066 AliTPCLaserTrack::LoadTracks();
1067 tree->SetBranchAddress("dY.",&vecdY);
1068 tree->SetBranchAddress("dZ.",&vecdZ);
1069 tree->SetBranchAddress("eY.",&veceY);
1070 tree->SetBranchAddress("eZ.",&veceZ);
1071 tree->SetBranchAddress("LTr.",<r);
1072 Int_t entries= tree->GetEntries();
1073 TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1074 Double_t bz=AliTrackerBase::GetBz();
1077 for (Int_t ientry=0; ientry<entries; ientry++){
1078 tree->GetEntry(ientry);
1079 if (!ltr->GetVecGX()){
1080 ltr->UpdatePoints();
1082 TVectorD * delta= (itype==0)? vecdY:vecdZ;
1083 TVectorD * err= (itype==0)? veceY:veceZ;
1085 for (Int_t irow=0; irow<159; irow++){
1086 Int_t nentries = 1000;
1087 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1088 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1090 Double_t phi =(*ltr->GetVecPhi())[irow];
1091 Double_t theta =ltr->GetTgl();
1092 Double_t mean=delta->GetMatrixArray()[irow];
1093 Double_t gx=0,gy=0,gz=0;
1094 Double_t snp = (*ltr->GetVecP2())[irow];
1095 Double_t rms = 0.1+err->GetMatrixArray()[irow];
1096 gx = (*ltr->GetVecGX())[irow];
1097 gy = (*ltr->GetVecGY())[irow];
1098 gz = (*ltr->GetVecGZ())[irow];
1099 Int_t bundle= ltr->GetBundle();
1102 // get delta R used in reconstruction
1103 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
1104 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1105 const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1106 Double_t xyz0[3]={gx,gy,gz};
1107 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1109 // old ExB correction
1111 if(recoParam&&recoParam->GetUseExBCorrection()) {
1112 Double_t xyz1[3]={gx,gy,gz};
1113 calib->GetExB()->Correct(xyz0,xyz1);
1114 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1117 if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1118 Float_t xyz1[3]={gx,gy,gz};
1119 Int_t sector=(gz>0)?0:18;
1120 correction->CorrectPoint(xyz1, sector);
1121 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1126 (*pcstream)<<"fit"<<
1127 "bz="<<bz<< // magnetic filed used
1128 "dtype="<<dtype<< // detector match type
1129 "ptype="<<itype<< // parameter type
1130 "theta="<<theta<< // theta
1131 "phi="<<phi<< // phi
1132 "snp="<<snp<< // snp
1133 "mean="<<mean<< // mean dist value
1134 "rms="<<rms<< // rms
1135 "gx="<<gx<< // global position
1136 "gy="<<gy<< // global position
1137 "gz="<<gz<< // global position
1138 "dRrec="<<dRrec<< // delta Radius in reconstruction
1139 "id="<<bundle<< //bundle
1140 "entries="<<nentries;// number of entries in bin
1143 Double_t ky = TMath::Tan(TMath::ASin(snp));
1144 Int_t ncorr = corrArray->GetEntries();
1145 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
1146 Double_t phi0 = TMath::ATan2(gy,gx);
1147 Double_t distortions[1000]={0};
1148 Double_t distortionsR[1000]={0};
1149 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1150 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1151 Float_t distPoint[3]={gx,gy,gz};
1152 Int_t sector= (gz>0)? 0:18;
1154 corr->DistortPoint(distPoint, sector);
1156 // Double_t value=distPoint[2]-gz;
1158 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1159 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1160 Double_t drphi= r0*(phi1-phi0);
1161 Double_t dr = r1-r0;
1162 distortions[icorr] = drphi-ky*dr;
1163 distortionsR[icorr] = dr;
1165 (*pcstream)<<"fit"<<
1166 Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
1167 Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
1169 (*pcstream)<<"fit"<<"\n";
1177 void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
1179 // make a distortion map out ou fthe residual histogram
1180 // Results are written to the debug streamer - pcstream
1182 // his0 - input (4D) residual histogram
1183 // pcstream - file to write the tree
1185 // marian.ivanov@cern.ch
1186 const Int_t kMinEntries=50;
1187 Int_t nbins1=his0->GetAxis(1)->GetNbins();
1188 Int_t first1=his0->GetAxis(1)->GetFirst();
1189 Int_t last1 =his0->GetAxis(1)->GetLast();
1191 Double_t bz=AliTrackerBase::GetBz();
1192 Int_t idim[4]={0,1,2,3};
1193 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
1195 his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1196 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1197 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
1198 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
1199 Int_t first3 = his1->GetAxis(3)->GetFirst();
1200 Int_t last3 = his1->GetAxis(3)->GetLast();
1203 for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){ // axis 3 - local angle
1204 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1205 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1207 his1->GetAxis(3)->SetRangeUser(-1,1);
1210 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
1211 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1212 Int_t first2 = his3->GetAxis(2)->GetFirst();
1213 Int_t last2 = his3->GetAxis(2)->GetLast();
1215 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1216 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1217 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1218 TH1 * hisDelta = his3->Projection(0);
1220 Double_t entries = hisDelta->GetEntries();
1221 Double_t mean=0, rms=0;
1222 if (entries>kMinEntries){
1223 mean = hisDelta->GetMean();
1224 rms = hisDelta->GetRMS();
1226 (*pcstream)<<hname<<
1232 "entries="<<entries<<
1237 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1249 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1251 // Store object in the OCDB
1252 // By default the object is stored in the current directory
1253 // default comment consit of user name and the date
1255 TString ocdbStorage="";
1256 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1257 AliCDBMetaData *metaData= new AliCDBMetaData();
1258 metaData->SetObjectClassName("AliTPCCorrection");
1259 metaData->SetResponsible("Marian Ivanov");
1260 metaData->SetBeamPeriod(1);
1261 metaData->SetAliRootVersion("05-25-01"); //root version
1262 TString userName=gSystem->GetFromPipe("echo $USER");
1263 TString date=gSystem->GetFromPipe("date");
1265 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1266 if (comment) metaData->SetComment(comment);
1268 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1269 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1270 gStorage->Put(this, (*id1), metaData);
1274 void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
1276 AliVertexerTracks *vertexer = new AliVertexerTracks(5);// 5kGaus
1278 TObjArray ATrk; // Original Track array of Aside
1279 TObjArray dATrk; // Distorted Track array of A side
1280 UShort_t *AId = new UShort_t[nTracks]; // A side Track ID
1283 UShort_t *CId = new UShort_t [nTracks];
1285 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1286 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
1287 fpt.SetParameters(7.24,0.120);
1289 for(Int_t nt=0; nt<nTracks; nt++){
1290 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
1291 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
1292 Double_t pt = fpt.GetRandom();// momentum for f1
1294 if(gRandom->Rndm() < 0.5){
1300 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1302 pxyz[0]=pt*TMath::Cos(phi);
1303 pxyz[1]=pt*TMath::Sin(phi);
1304 pxyz[2]=pt*TMath::Tan(theta);
1305 Double_t cv[21]={0};
1306 AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
1310 AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir, NULL);
1312 if (pcstream) (*pcstream)<<"track"<<
1318 if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
1322 Int_t nn=ATrk.GetEntriesFast();
1325 }else if(( eta<-0.07 )&&( eta>-etaCuts )){
1329 Int_t nn=CTrk.GetEntriesFast();
1334 }// end of track loop
1336 vertexer->SetTPCMode();
1337 vertexer->SetConstraintOff();
1339 aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dATrk,AId));
1340 avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&ATrk,AId));
1341 cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dCTrk,CId));
1342 cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&CTrk,CId));
1343 if (pcstream) (*pcstream)<<"vertex"<<
1344 "x="<<orgVertex[0]<<
1345 "y="<<orgVertex[1]<<
1346 "z="<<orgVertex[2]<<
1347 "av.="<<&aV<< // distorted vertex A side
1348 "cv.="<<&cV<< // distroted vertex C side
1349 "avO.="<<&avOrg<< // original vertex A side
1356 void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
1358 // make correction available for visualization using
1359 // TFormula, TFX and TTree::Draw
1360 // important in order to check corrections and also compute dervied variables
1361 // e.g correction partial derivatives
1363 // NOTE - class is not owner of correction
1365 if (!fgVisualCorrection) fgVisualCorrection=new TObjArray;
1366 fgVisualCorrection->AddAt(corr, position);
1371 Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t localX, Double_t kZ, Int_t axisType, Int_t corrType){
1373 // calculate the correction at given position - check the geffCorr
1375 if (!fgVisualCorrection) return 0;
1376 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1377 if (!corr) return 0;
1378 Double_t phi=sector*TMath::Pi()/9.;
1379 Double_t gx = localX*TMath::Cos(phi);
1380 Double_t gy = localX*TMath::Sin(phi);
1381 Double_t gz = localX*kZ;
1382 Int_t nsector=(gz>0) ? 0:18;
1386 Float_t distPoint[3]={gx,gy,gz};
1387 corr->DistortPoint(distPoint, nsector);
1388 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1389 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1390 Double_t phi0=TMath::ATan2(gy,gx);
1391 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1392 if (axisType==0) return r1-r0;
1393 if (axisType==1) return (phi1-phi0)*r0;
1394 if (axisType==2) return distPoint[2]-gz;
1398 Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
1400 // return correction at given x,y,z
1402 if (!fgVisualCorrection) return 0;
1403 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1404 if (!corr) return 0;
1405 Double_t phi0= TMath::ATan2(gy,gx);
1406 Int_t nsector=(gz>0) ? 0:18;
1407 Float_t distPoint[3]={gx,gy,gz};
1408 corr->DistortPoint(distPoint, nsector);
1409 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1410 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1411 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1412 if (axisType==0) return r1-r0;
1413 if (axisType==1) return (phi1-phi0)*r0;
1414 if (axisType==2) return distPoint[2]-gz;