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>
54 #include "AliTPCParamSR.h"
56 #include "AliTPCCorrection.h"
59 #include "AliExternalTrackParam.h"
60 #include "AliTrackPointArray.h"
61 #include "TDatabasePDG.h"
62 #include "AliTrackerBase.h"
63 #include "AliTPCROC.h"
64 #include "THnSparse.h"
66 #include "AliTPCLaserTrack.h"
67 #include "AliESDVertex.h"
68 #include "AliVertexerTracks.h"
69 #include "TDatabasePDG.h"
73 #include "TDatabasePDG.h"
75 #include "AliTPCTransform.h"
76 #include "AliTPCcalibDB.h"
77 #include "AliTPCExB.h"
79 #include "AliTPCRecoParam.h"
82 ClassImp(AliTPCCorrection)
85 TObjArray *AliTPCCorrection::fgVisualCorrection=0;
86 // instance of correction for visualization
89 // FIXME: the following values should come from the database
90 const Double_t AliTPCCorrection::fgkTPCZ0 = 249.7; // nominal gating grid position
91 const Double_t AliTPCCorrection::fgkIFCRadius= 83.5; // radius which renders the "18 rod manifold" best -> compare calc. of Jim Thomas
92 // compare gkIFCRadius= 83.05: Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
93 const Double_t AliTPCCorrection::fgkOFCRadius= 254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
94 const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
95 const Double_t AliTPCCorrection::fgkCathodeV = -100000.0; // Cathode Voltage (volts)
96 const Double_t AliTPCCorrection::fgkGG = -70.0; // Gating Grid voltage (volts)
98 const Double_t AliTPCCorrection::fgkdvdE = 0.0024; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
100 const Double_t AliTPCCorrection::fgkEM = -1.602176487e-19/9.10938215e-31; // charge/mass in [C/kg]
101 const Double_t AliTPCCorrection::fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)]
103 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
104 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
105 83.06, 83.5, 84.0, 84.5, 85.0, 85.5, 86.0, 86.5,
106 87.0, 87.5, 88.0, 88.5, 89.0, 89.5, 90.0, 90.5, 91.0, 92.0,
107 93.0, 94.0, 95.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, 241.0, 242.0, 243.0, 244.0,
115 245.0, 245.5, 246.0, 246.5, 247.0, 247.5, 248.0, 248.5, 249.0, 249.5,
116 250.0, 250.5, 251.0, 251.5, 252.0, 252.5, 253.0, 253.5, 254.0, 254.5
119 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = {
120 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
121 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
122 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
123 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
124 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
125 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
126 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
127 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
128 -100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0,
129 -80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0,
130 -60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0,
131 -40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0,
132 -20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0,
133 -1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0,
134 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
135 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0,
136 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0,
137 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0,
138 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0,
139 102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
140 122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
141 142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
142 162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
143 182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
144 202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0,
145 222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0,
146 242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ;
150 AliTPCCorrection::AliTPCCorrection()
151 : TNamed("correction_unity","unity"),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
154 // default constructor
156 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
158 // Initialization of interpolation points
159 for (Int_t i = 0; i<kNPhi; i++) {
160 fgkPhiList[i] = TMath::TwoPi()*i/(kNPhi-1);
165 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
166 : TNamed(name,title),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
169 // default constructor, that set the name and title
171 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
173 // Initialization of interpolation points
174 for (Int_t i = 0; i<kNPhi; i++) {
175 fgkPhiList[i] = TMath::TwoPi()*i/(kNPhi-1);
180 AliTPCCorrection::~AliTPCCorrection() {
182 // virtual destructor
186 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
188 // Corrects the initial coordinates x (cartesian coordinates)
189 // according to the given effect (inherited classes)
190 // roc represents the TPC read out chamber (offline numbering convention)
193 GetCorrection(x,roc,dx);
194 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
197 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
199 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
200 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
201 // roc represents the TPC read out chamber (offline numbering convention)
204 GetCorrection(x,roc,dx);
205 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
208 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
210 // Distorts the initial coordinates x (cartesian coordinates)
211 // according to the given effect (inherited classes)
212 // roc represents the TPC read out chamber (offline numbering convention)
215 GetDistortion(x,roc,dx);
216 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
219 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
221 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
222 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
223 // roc represents the TPC read out chamber (offline numbering convention)
226 GetDistortion(x,roc,dx);
227 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
230 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
232 // This function delivers the correction values dx in respect to the inital coordinates x
233 // roc represents the TPC read out chamber (offline numbering convention)
234 // Note: The dx is overwritten by the inherited effectice class ...
236 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
239 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
241 // This function delivers the distortion values dx in respect to the inital coordinates x
242 // roc represents the TPC read out chamber (offline numbering convention)
244 GetCorrection(x,roc,dx);
245 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
248 void AliTPCCorrection::Init() {
250 // Initialization funtion (not used at the moment)
254 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
260 void AliTPCCorrection::Print(Option_t* /*option*/) const {
262 // Print function to check which correction classes are used
263 // option=="d" prints details regarding the setted magnitude
264 // option=="a" prints the C0 and C1 coefficents for calibration purposes
266 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
269 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
271 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
272 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
277 //SetOmegaTauT1T2(omegaTau, t1, t2);
280 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
282 // Simple plot functionality.
283 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
284 // in respect to position z within the XY plane.
285 // The histogramm has nx times ny entries.
287 AliTPCParam* tpcparam = new AliTPCParamSR;
289 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
290 nx,-250.,250.,ny,-250.,250.);
293 Int_t roc=z>0.?0:18; // FIXME
294 for (Int_t iy=1;iy<=ny;++iy) {
295 x[1]=h->GetYaxis()->GetBinCenter(iy);
296 for (Int_t ix=1;ix<=nx;++ix) {
297 x[0]=h->GetXaxis()->GetBinCenter(ix);
298 GetCorrection(x,roc,dx);
299 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
300 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
301 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
302 h->SetBinContent(ix,iy,r1-r0);
305 h->SetBinContent(ix,iy,0.);
312 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
314 // Simple plot functionality.
315 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
316 // in respect to position z within the XY plane.
317 // The histogramm has nx times ny entries.
320 AliTPCParam* tpcparam = new AliTPCParamSR;
322 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
323 nx,-250.,250.,ny,-250.,250.);
326 Int_t roc=z>0.?0:18; // FIXME
327 for (Int_t iy=1;iy<=ny;++iy) {
328 x[1]=h->GetYaxis()->GetBinCenter(iy);
329 for (Int_t ix=1;ix<=nx;++ix) {
330 x[0]=h->GetXaxis()->GetBinCenter(ix);
331 GetCorrection(x,roc,dx);
332 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
333 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
334 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
335 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
337 Float_t dphi=phi1-phi0;
338 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
339 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
341 h->SetBinContent(ix,iy,r0*dphi);
344 h->SetBinContent(ix,iy,0.);
351 TH2F* AliTPCCorrection::CreateHistoDZinXY(Float_t z,Int_t nx,Int_t ny) {
353 // Simple plot functionality.
354 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
355 // in respect to position z within the XY plane.
356 // The histogramm has nx times ny entries.
359 AliTPCParam* tpcparam = new AliTPCParamSR;
361 TH2F *h=CreateTH2F("dz_xy",GetTitle(),"x [cm]","y [cm]","dz [cm]",
362 nx,-250.,250.,ny,-250.,250.);
365 Int_t roc=z>0.?0:18; // FIXME
366 for (Int_t iy=1;iy<=ny;++iy) {
367 x[1]=h->GetYaxis()->GetBinCenter(iy);
368 for (Int_t ix=1;ix<=nx;++ix) {
369 x[0]=h->GetXaxis()->GetBinCenter(ix);
370 GetCorrection(x,roc,dx);
371 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
372 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
373 h->SetBinContent(ix,iy,dx[2]);
376 h->SetBinContent(ix,iy,0.);
383 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
385 // Simple plot functionality.
386 // Returns a 2d hisogram which represents the corrections in r direction (dr)
387 // in respect to angle phi within the ZR plane.
388 // The histogramm has nx times ny entries.
390 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
391 nz,-250.,250.,nr,85.,250.);
393 for (Int_t ir=1;ir<=nr;++ir) {
394 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
395 x[0]=radius*TMath::Cos(phi);
396 x[1]=radius*TMath::Sin(phi);
397 for (Int_t iz=1;iz<=nz;++iz) {
398 x[2]=h->GetXaxis()->GetBinCenter(iz);
399 Int_t roc=x[2]>0.?0:18; // FIXME
400 GetCorrection(x,roc,dx);
401 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
402 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
403 h->SetBinContent(iz,ir,r1-r0);
410 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
412 // Simple plot functionality.
413 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
414 // in respect to angle phi within the ZR plane.
415 // The histogramm has nx times ny entries.
417 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
418 nz,-250.,250.,nr,85.,250.);
420 for (Int_t iz=1;iz<=nz;++iz) {
421 x[2]=h->GetXaxis()->GetBinCenter(iz);
422 Int_t roc=x[2]>0.?0:18; // FIXME
423 for (Int_t ir=1;ir<=nr;++ir) {
424 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
425 x[0]=radius*TMath::Cos(phi);
426 x[1]=radius*TMath::Sin(phi);
427 GetCorrection(x,roc,dx);
428 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
429 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
430 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
432 Float_t dphi=phi1-phi0;
433 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
434 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
436 h->SetBinContent(iz,ir,r0*dphi);
442 TH2F* AliTPCCorrection::CreateHistoDZinZR(Float_t phi,Int_t nz,Int_t nr) {
444 // Simple plot functionality.
445 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
446 // in respect to angle phi within the ZR plane.
447 // The histogramm has nx times ny entries.
449 TH2F *h=CreateTH2F("dz_zr",GetTitle(),"z [cm]","r [cm]","dz [cm]",
450 nz,-250.,250.,nr,85.,250.);
452 for (Int_t ir=1;ir<=nr;++ir) {
453 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
454 x[0]=radius*TMath::Cos(phi);
455 x[1]=radius*TMath::Sin(phi);
456 for (Int_t iz=1;iz<=nz;++iz) {
457 x[2]=h->GetXaxis()->GetBinCenter(iz);
458 Int_t roc=x[2]>0.?0:18; // FIXME
459 GetCorrection(x,roc,dx);
460 h->SetBinContent(iz,ir,dx[2]);
468 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
469 const char *xlabel,const char *ylabel,const char *zlabel,
470 Int_t nbinsx,Double_t xlow,Double_t xup,
471 Int_t nbinsy,Double_t ylow,Double_t yup) {
473 // Helper function to create a 2d histogramm of given size
479 while (gDirectory->FindObject(hname.Data())) {
486 TH2F *h=new TH2F(hname.Data(),title,
489 h->GetXaxis()->SetTitle(xlabel);
490 h->GetYaxis()->SetTitle(ylabel);
491 h->GetZaxis()->SetTitle(zlabel);
496 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
498 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
499 const Double_t er[kNZ][kNR], Double_t &erValue ) {
501 // Interpolate table - 2D interpolation
503 Double_t saveEr[5] = {0,0,0,0,0};
505 Search( kNZ, fgkZList, z, fJLow ) ;
506 Search( kNR, fgkRList, r, fKLow ) ;
507 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
508 if ( fKLow < 0 ) fKLow = 0 ;
509 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
510 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
512 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
513 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
515 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
519 void AliTPCCorrection::Interpolate3DEdistortion( const Int_t order, const Double_t r, const Float_t phi, const Double_t z,
520 const Double_t er[kNZ][kNPhi][kNR], const Double_t ephi[kNZ][kNPhi][kNR], const Double_t ez[kNZ][kNPhi][kNR],
521 Double_t &erValue, Double_t &ephiValue, Double_t &ezValue) {
523 // Interpolate table - 3D interpolation
526 Double_t saveEr[5]= {0,0,0,0,0};
527 Double_t savedEr[5]= {0,0,0,0,0} ;
529 Double_t saveEphi[5]= {0,0,0,0,0};
530 Double_t savedEphi[5]= {0,0,0,0,0} ;
532 Double_t saveEz[5]= {0,0,0,0,0};
533 Double_t savedEz[5]= {0,0,0,0,0} ;
535 Search( kNZ, fgkZList, z, fILow ) ;
536 Search( kNPhi, fgkPhiList, z, fJLow ) ;
537 Search( kNR, fgkRList, r, fKLow ) ;
539 if ( fILow < 0 ) fILow = 0 ; // check if out of range
540 if ( fJLow < 0 ) fJLow = 0 ;
541 if ( fKLow < 0 ) fKLow = 0 ;
543 if ( fILow + order >= kNZ - 1 ) fILow = kNZ - 1 - order ;
544 if ( fJLow + order >= kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
545 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
547 for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
548 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
549 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r ) ;
550 saveEphi[j-fJLow] = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
551 saveEz[j-fJLow] = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r ) ;
553 savedEr[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi ) ;
554 savedEphi[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
555 savedEz[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi ) ;
557 erValue = Interpolate( &fgkZList[fILow], savedEr, order, z ) ;
558 ephiValue = Interpolate( &fgkZList[fILow], savedEphi, order, z ) ;
559 ezValue = Interpolate( &fgkZList[fILow], savedEz, order, z ) ;
563 Double_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
564 const Int_t nx, const Int_t ny, const Double_t xv[], const Double_t yv[],
565 const TMatrixD &array ) {
567 // Interpolate table (TMatrix format) - 2D interpolation
570 static Int_t jlow = 0, klow = 0 ;
571 Double_t saveArray[5] = {0,0,0,0,0} ;
573 Search( nx, xv, x, jlow ) ;
574 Search( ny, yv, y, klow ) ;
575 if ( jlow < 0 ) jlow = 0 ; // check if out of range
576 if ( klow < 0 ) klow = 0 ;
577 if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
578 if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
580 for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
582 Double_t *ajkl = &((TMatrixD&)array)(j,klow);
583 saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
586 return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
590 Double_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x, const Double_t y, const Double_t z,
591 const Int_t nx, const Int_t ny, const Int_t nz,
592 const Double_t xv[], const Double_t yv[], const Double_t zv[],
593 TMatrixD **arrayofArrays ) {
595 // Interpolate table (TMatrix format) - 3D interpolation
598 static Int_t ilow = 0, jlow = 0, klow = 0 ;
599 Double_t saveArray[5]= {0,0,0,0,0};
600 Double_t savedArray[5]= {0,0,0,0,0} ;
602 Search( nx, xv, x, ilow ) ;
603 Search( ny, yv, y, jlow ) ;
604 Search( nz, zv, z, klow ) ;
606 if ( ilow < 0 ) ilow = 0 ; // check if out of range
607 if ( jlow < 0 ) jlow = 0 ;
608 if ( klow < 0 ) klow = 0 ;
610 if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
611 if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
612 if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
614 for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
616 TMatrixD &table = *arrayofArrays[k] ;
617 for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
619 saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
621 savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
623 return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
629 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
630 const Int_t order, const Double_t x ) {
632 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
636 if ( order == 2 ) { // Quadratic Interpolation = 2
637 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
638 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
639 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
640 } else { // Linear Interpolation = 1
641 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
649 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
651 // Search an ordered table by starting at the most recently used point
654 Long_t middle, high ;
655 Int_t ascend = 0, increment = 1 ;
657 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
659 if ( low < 0 || low > n-1 ) {
660 low = -1 ; high = n ;
661 } else { // Ordered Search phase
662 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
663 if ( low == n-1 ) return ;
665 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
668 high = low + increment ;
669 if ( high > n-1 ) { high = n ; break ; }
672 if ( low == 0 ) { low = -1 ; return ; }
674 while ( (Int_t)( x < xArray[low] ) == ascend ) {
677 if ( increment >= high ) { low = -1 ; break ; }
678 else low = high - increment ;
683 while ( (high-low) != 1 ) { // Binary Search Phase
684 middle = ( high + low ) / 2 ;
685 if ( (Int_t)( x >= xArray[middle] ) == ascend )
691 if ( x == xArray[n-1] ) low = n-2 ;
692 if ( x == xArray[0] ) low = 0 ;
696 void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
697 TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
698 const Int_t rows, const Int_t columns, const Int_t iterations,
699 const Bool_t rocDisplacement ) {
701 // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
703 // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
704 // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
705 // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
706 // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
707 // you wish to solve Laplaces equation however it should not contain random numbers or you will get
708 // random numbers back as a solution.
709 // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
710 // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
711 // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
712 // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
713 // number of points and solves the problem again. This happens several times until the maximum number
714 // of points has been included in the array.
716 // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
717 // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
719 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
721 // Original code by Jim Thomas (STAR TPC Collaboration)
724 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
726 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
727 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
728 const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
730 TMatrixD arrayEr(rows,columns) ;
731 TMatrixD arrayEz(rows,columns) ;
733 //Check that number of rows and columns is suitable for a binary expansion
735 if ( !IsPowerOfTwo(rows-1) ) {
736 AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
739 if ( !IsPowerOfTwo(columns-1) ) {
740 AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
744 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
745 // Allow for different size grid spacing in R and Z directions
746 // Use a binary expansion of the size of the matrix to speed up the solution of the problem
748 Int_t iOne = (rows-1)/4 ;
749 Int_t jOne = (columns-1)/4 ;
750 // Solve for N in 2**N, add one.
751 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
753 for ( Int_t count = 0 ; count < loops ; count++ ) {
754 // Loop while the matrix expands & the resolution increases.
756 Float_t tempGridSizeR = gridSizeR * iOne ;
757 Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
758 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
760 // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
761 std::vector<float> coef1(rows) ;
762 std::vector<float> coef2(rows) ;
764 for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
765 Float_t radius = fgkIFCRadius + i*gridSizeR ;
766 coef1[i] = 1.0 + tempGridSizeR/(2*radius);
767 coef2[i] = 1.0 - tempGridSizeR/(2*radius);
770 TMatrixD sumChargeDensity(rows,columns) ;
772 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
773 Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
774 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
775 if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
777 // Add up all enclosed charge density contributions within 1/2 unit in all directions
778 Float_t weight = 0.0 ;
780 sumChargeDensity(i,j) = 0.0 ;
781 for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
782 for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
783 if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
786 // Note that this is cylindrical geometry
787 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
788 sum += weight*radius ;
791 sumChargeDensity(i,j) /= sum ;
793 sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
797 for ( Int_t k = 1 ; k <= iterations; k++ ) {
798 // Solve Poisson's Equation
799 // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
800 // as interations increase.
801 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
802 Float_t overRelaxM1 = overRelax - 1.0 ;
803 Float_t overRelaxtempFourth, overRelaxcoef5 ;
804 overRelaxtempFourth = overRelax * tempFourth ;
805 overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
807 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
808 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
810 arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
811 + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
812 - overRelaxcoef5 * arrayV(i,j)
813 + coef1[i] * arrayV(i+iOne,j)
814 + sumChargeDensity(i,j)
815 ) * overRelaxtempFourth;
819 if ( k == iterations ) {
820 // After full solution is achieved, copy low resolution solution into higher res array
821 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
822 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
825 arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
826 if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
829 arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
830 if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
832 if ( iOne > 1 && jOne > 1 ) {
833 arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
834 if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
835 if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
836 // Note that this leaves a point at the upper left and lower right corners uninitialized.
837 // -> Not a big deal.
846 iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
847 jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
849 sumChargeDensity.Clear();
852 // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
853 for ( Int_t j = 0 ; j < columns ; j++ ) {
854 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
855 arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
856 arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
859 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
860 for ( Int_t i = 0 ; i < rows ; i++) {
861 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
862 arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
863 arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
866 for ( Int_t i = 0 ; i < rows ; i++) {
867 // Note: go back and compare to old version of this code. See notes below.
868 // JT Test ... attempt to divide by real Ez not Ez to first order
869 for ( Int_t j = 0 ; j < columns ; j++ ) {
870 arrayEz(i,j) += ezField;
871 // This adds back the overall Z gradient of the field (main E field component)
873 // Warning: (-=) assumes you are using an error potetial without the overall Field included
876 // Integrate Er/Ez from Z to zero
877 for ( Int_t j = 0 ; j < columns ; j++ ) {
878 for ( Int_t i = 0 ; i < rows ; i++ ) {
880 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
881 arrayErOverEz(i,j) = 0.0 ;
882 arrayDeltaEz(i,j) = 0.0 ;
884 for ( Int_t k = j ; k < columns ; k++ ) {
885 arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
886 arrayDeltaEz(i,j) += index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
887 if ( index != 4 ) index = 4; else index = 2 ;
890 arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
891 arrayDeltaEz(i,j) -= (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
894 arrayErOverEz(i,j) += (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
895 -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
896 arrayDeltaEz(i,j) += (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
897 -2.5*(arrayEz(i,columns-1)-ezField));
899 if ( j == columns-2 ) {
900 arrayErOverEz(i,j) = (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
901 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
902 arrayDeltaEz(i,j) = (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
903 +1.5*(arrayEz(i,columns-1)-ezField) ) ;
905 if ( j == columns-1 ) {
906 arrayErOverEz(i,j) = 0.0 ;
907 arrayDeltaEz(i,j) = 0.0 ;
912 // calculate z distortion from the integrated Delta Ez residuals
913 // and include the aquivalence (Volt to cm) of the ROC shift !!
915 for ( Int_t j = 0 ; j < columns ; j++ ) {
916 for ( Int_t i = 0 ; i < rows ; i++ ) {
918 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
919 arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
921 // ROC Potential in cm aquivalent
922 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
923 if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift; // add the ROC misaligment
933 void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
934 TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
935 const Int_t rows, const Int_t columns, const Int_t phislices,
936 const Float_t deltaphi, const Int_t iterations, const Int_t symmetry,
937 Bool_t rocDisplacement ) {
939 // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
941 // NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
942 // The number of rows and COLUMNS can be different.
945 // COLUMNS == 2**N + 1
946 // PHISLICES == Arbitrary but greater than 3
948 // DeltaPhi in Radians
950 // SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
951 // = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
953 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
955 const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
957 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
958 const Float_t gridSizePhi = deltaphi ;
959 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
960 const Float_t ratioPhi = gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
961 const Float_t ratioZ = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
963 TMatrixD arrayE(rows,columns) ;
965 // Check that the number of rows and columns is suitable for a binary expansion
966 if ( !IsPowerOfTwo((rows-1)) ) {
967 AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
969 if ( !IsPowerOfTwo((columns-1)) ) {
970 AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
972 if ( phislices <= 3 ) {
973 AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
975 if ( phislices > 1000 ) {
976 AliError("Poisson3D phislices > 1000 is not allowed (nor wise) ");
979 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
980 // Allow for different size grid spacing in R and Z directions
981 // Use a binary expansion of the matrix to speed up the solution of the problem
983 Int_t loops, mplus, mminus, signplus, signminus ;
984 Int_t ione = (rows-1)/4 ;
985 Int_t jone = (columns-1)/4 ;
986 loops = TMath::Max(ione, jone) ; // Calculate the number of loops for the binary expansion
987 loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ; // Solve for N in 2**N
989 TMatrixD* arrayofSumChargeDensities[1000] ; // Create temporary arrays to store low resolution charge arrays
991 for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
993 for ( Int_t count = 0 ; count < loops ; count++ ) { // START the master loop and do the binary expansion
995 Float_t tempgridSizeR = gridSizeR * ione ;
996 Float_t tempratioPhi = ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
997 Float_t tempratioZ = ratioZ * ione * ione / ( jone * jone ) ;
999 std::vector<float> coef1(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1000 std::vector<float> coef2(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1001 std::vector<float> coef3(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1002 std::vector<float> coef4(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1004 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1005 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1006 coef1[i] = 1.0 + tempgridSizeR/(2*radius);
1007 coef2[i] = 1.0 - tempgridSizeR/(2*radius);
1008 coef3[i] = tempratioPhi/(radius*radius);
1009 coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
1012 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1013 TMatrixD &chargeDensity = *arrayofChargeDensities[m] ;
1014 TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
1015 for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
1016 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1017 for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
1018 if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1019 else { // Add up all enclosed charge density contributions within 1/2 unit in all directions
1020 Float_t weight = 0.0 ;
1022 sumChargeDensity(i,j) = 0.0 ;
1023 for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
1024 for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
1025 if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
1028 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1029 sum += weight*radius ;
1032 sumChargeDensity(i,j) /= sum ;
1034 sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
1039 for ( Int_t k = 1 ; k <= iterations; k++ ) {
1041 // over-relaxation index, >= 1 but < 2
1042 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1043 Float_t overRelaxM1 = overRelax - 1.0 ;
1045 std::vector<float> overRelaxcoef4(rows) ; // Do this the standard C++ way to avoid gcc extensions
1046 std::vector<float> overRelaxcoef5(rows) ; // Do this the standard C++ way to avoid gcc extensions
1048 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1049 overRelaxcoef4[i] = overRelax * coef4[i] ;
1050 overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
1053 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1055 mplus = m + 1; signplus = 1 ;
1056 mminus = m - 1 ; signminus = 1 ;
1057 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1058 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1059 if ( mminus < 0 ) mminus = 1 ;
1061 else if (symmetry==-1) { // Anti-symmetry in phi
1062 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1063 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1065 else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1066 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1067 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1069 TMatrixD& arrayV = *arrayofArrayV[m] ;
1070 TMatrixD& arrayVP = *arrayofArrayV[mplus] ;
1071 TMatrixD& arrayVM = *arrayofArrayV[mminus] ;
1072 TMatrixD& sumChargeDensity = *arrayofSumChargeDensities[m] ;
1074 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1075 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1077 arrayV(i,j) = ( coef2[i] * arrayV(i-ione,j)
1078 + tempratioZ * ( arrayV(i,j-jone) + arrayV(i,j+jone) )
1079 - overRelaxcoef5[i] * arrayV(i,j)
1080 + coef1[i] * arrayV(i+ione,j)
1081 + coef3[i] * ( signplus*arrayVP(i,j) + signminus*arrayVM(i,j) )
1082 + sumChargeDensity(i,j)
1083 ) * overRelaxcoef4[i] ;
1084 // Note: over-relax the solution at each step. This speeds up the convergance.
1089 if ( k == iterations ) { // After full solution is achieved, copy low resolution solution into higher res array
1090 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1091 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1094 arrayV(i+ione/2,j) = ( arrayV(i+ione,j) + arrayV(i,j) ) / 2 ;
1095 if ( i == ione ) arrayV(i-ione/2,j) = ( arrayV(0,j) + arrayV(ione,j) ) / 2 ;
1098 arrayV(i,j+jone/2) = ( arrayV(i,j+jone) + arrayV(i,j) ) / 2 ;
1099 if ( j == jone ) arrayV(i,j-jone/2) = ( arrayV(i,0) + arrayV(i,jone) ) / 2 ;
1101 if ( ione > 1 && jone > 1 ) {
1102 arrayV(i+ione/2,j+jone/2) = ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
1103 if ( i == ione ) arrayV(i-ione/2,j-jone/2) = ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
1104 if ( j == jone ) arrayV(i-ione/2,j-jone/2) = ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
1105 // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
1114 ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
1115 jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
1119 //Differentiate V(r) and solve for E(r) using special equations for the first and last row
1120 //Integrate E(r)/E(z) from point of origin to pad plane
1122 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1123 TMatrixD& arrayV = *arrayofArrayV[m] ;
1124 TMatrixD& eroverEz = *arrayofEroverEz[m] ;
1126 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1128 // Differentiate in R
1129 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1130 arrayE(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1131 arrayE(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1133 for ( Int_t i = 0 ; i < rows ; i++ ) {
1134 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1135 eroverEz(i,j) = 0.0 ;
1136 for ( Int_t k = j ; k < columns ; k++ ) {
1138 eroverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1139 if ( index != 4 ) index = 4; else index = 2 ;
1141 if ( index == 4 ) eroverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1142 if ( index == 2 ) eroverEz(i,j) +=
1143 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1144 if ( j == columns-2 ) eroverEz(i,j) =
1145 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1146 if ( j == columns-1 ) eroverEz(i,j) = 0.0 ;
1149 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1150 // eroverEz.Draw("surf") ; } // JT test
1153 //Differentiate V(r) and solve for E(phi)
1154 //Integrate E(phi)/E(z) from point of origin to pad plane
1156 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1158 mplus = m + 1; signplus = 1 ;
1159 mminus = m - 1 ; signminus = 1 ;
1160 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1161 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1162 if ( mminus < 0 ) mminus = 1 ;
1164 else if (symmetry==-1) { // Anti-symmetry in phi
1165 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1166 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1168 else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
1169 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1170 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1172 TMatrixD &arrayVP = *arrayofArrayV[mplus] ;
1173 TMatrixD &arrayVM = *arrayofArrayV[mminus] ;
1174 TMatrixD &ePhioverEz = *arrayofEPhioverEz[m] ;
1175 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1176 // Differentiate in Phi
1177 for ( Int_t i = 0 ; i < rows ; i++ ) {
1178 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1179 arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
1182 for ( Int_t i = 0 ; i < rows ; i++ ) {
1183 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1184 ePhioverEz(i,j) = 0.0 ;
1185 for ( Int_t k = j ; k < columns ; k++ ) {
1187 ePhioverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1188 if ( index != 4 ) index = 4; else index = 2 ;
1190 if ( index == 4 ) ePhioverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1191 if ( index == 2 ) ePhioverEz(i,j) +=
1192 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1193 if ( j == columns-2 ) ePhioverEz(i,j) =
1194 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1195 if ( j == columns-1 ) ePhioverEz(i,j) = 0.0 ;
1198 // if ( m == 5 ) { TCanvas* c2 = new TCanvas("arrayE","arrayE",50,50,840,600) ; c2 -> cd() ;
1199 // arrayE.Draw("surf") ; } // JT test
1203 // Differentiate V(r) and solve for E(z) using special equations for the first and last row
1204 // Integrate (E(z)-Ezstd) from point of origin to pad plane
1206 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1207 TMatrixD& arrayV = *arrayofArrayV[m] ;
1208 TMatrixD& deltaEz = *arrayofDeltaEz[m] ;
1210 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1211 for ( Int_t i = 0 ; i < rows ; i++) {
1212 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1213 arrayE(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1214 arrayE(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1217 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1219 for ( Int_t i = 0 ; i < rows ; i++ ) {
1220 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1221 deltaEz(i,j) = 0.0 ;
1222 for ( Int_t k = j ; k < columns ; k++ ) {
1223 deltaEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k) ;
1224 if ( index != 4 ) index = 4; else index = 2 ;
1226 if ( index == 4 ) deltaEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1) ;
1227 if ( index == 2 ) deltaEz(i,j) +=
1228 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
1229 if ( j == columns-2 ) deltaEz(i,j) =
1230 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
1231 if ( j == columns-1 ) deltaEz(i,j) = 0.0 ;
1234 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1235 // eroverEz.Draw("surf") ; } // JT test
1237 // calculate z distortion from the integrated Delta Ez residuals
1238 // and include the aquivalence (Volt to cm) of the ROC shift !!
1240 for ( Int_t j = 0 ; j < columns ; j++ ) {
1241 for ( Int_t i = 0 ; i < rows ; i++ ) {
1243 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1244 deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
1246 // ROC Potential in cm aquivalent
1247 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
1248 if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift; // add the ROC misaligment
1253 } // end loop over phi
1257 for ( Int_t k = 0 ; k < phislices ; k++ )
1259 arrayofSumChargeDensities[k]->Delete() ;
1268 Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1270 // Helperfunction: Check if integer is a power of 2
1273 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
1274 if ( j == 1 ) return(1) ; // True
1275 return(0) ; // False
1279 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
1281 // Fit the track parameters - without and with distortion
1282 // 1. Space points in the TPC are simulated along the trajectory
1283 // 2. Space points distorted
1284 // 3. Fits the non distorted and distroted track to the reference plane at refX
1285 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
1287 // trackIn - input track parameters
1288 // refX - reference X to fit the track
1289 // dir - direction - out=1 or in=-1
1290 // pcstream - debug streamer to check the results
1292 // see AliExternalTrackParam.h documentation:
1293 // track1.fP[0] - local y (rphi)
1295 // track1.fP[2] - sinus of local inclination angle
1296 // track1.fP[3] - tangent of deep angle
1297 // track1.fP[4] - 1/pt
1299 AliTPCROC * roc = AliTPCROC::Instance();
1300 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
1301 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
1302 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1303 const Double_t kMaxSnp = 0.85;
1304 const Double_t kSigmaY=0.1;
1305 const Double_t kSigmaZ=0.1;
1306 const Double_t kMaxR=500;
1307 const Double_t kMaxZ=500;
1308 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1312 AliExternalTrackParam track(trackIn); //
1314 AliTrackPointArray pointArray0(npoints0);
1315 AliTrackPointArray pointArray1(npoints0);
1317 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
1319 // simulate the track
1321 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
1322 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
1323 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
1325 xyz[0]+=gRandom->Gaus(0,0.00005);
1326 xyz[1]+=gRandom->Gaus(0,0.00005);
1327 xyz[2]+=gRandom->Gaus(0,0.00005);
1328 if (TMath::Abs(track.GetZ())>kMaxZ) break;
1329 if (TMath::Abs(track.GetX())>kMaxR) break;
1330 AliTrackPoint pIn0; // space point
1332 Int_t sector= (xyz[2]>0)? 0:18;
1333 pointArray0.GetPoint(pIn0,npoints);
1334 pointArray1.GetPoint(pIn1,npoints);
1335 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
1336 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
1337 DistortPoint(distPoint, sector);
1338 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
1339 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
1341 track.Rotate(alpha);
1342 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
1343 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
1344 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
1345 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
1346 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
1347 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
1348 pointArray0.AddPoint(npoints, &pIn0);
1349 pointArray1.AddPoint(npoints, &pIn1);
1351 if (npoints>=npoints0) break;
1353 if (npoints<npoints0/2) return 0;
1357 AliExternalTrackParam *track0=0;
1358 AliExternalTrackParam *track1=0;
1359 AliTrackPoint point1,point2,point3;
1360 if (dir==1) { //make seed inner
1361 pointArray0.GetPoint(point1,1);
1362 pointArray0.GetPoint(point2,30);
1363 pointArray0.GetPoint(point3,60);
1365 if (dir==-1){ //make seed outer
1366 pointArray0.GetPoint(point1,npoints-60);
1367 pointArray0.GetPoint(point2,npoints-30);
1368 pointArray0.GetPoint(point3,npoints-1);
1370 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
1371 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
1373 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
1374 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
1378 pointArray0.GetPoint(pIn0,ipoint);
1379 pointArray1.GetPoint(pIn1,ipoint);
1380 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
1381 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
1383 if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1384 if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1385 if (TMath::Abs(track0->GetZ())>kMaxZ) break;
1386 if (TMath::Abs(track0->GetX())>kMaxR) break;
1387 if (TMath::Abs(track1->GetZ())>kMaxZ) break;
1388 if (TMath::Abs(track1->GetX())>kMaxR) break;
1390 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
1392 Double_t pointPos[2]={0,0};
1393 Double_t pointCov[3]={0,0,0};
1394 pointPos[0]=prot0.GetY();//local y
1395 pointPos[1]=prot0.GetZ();//local z
1396 pointCov[0]=prot0.GetCov()[3];//simay^2
1397 pointCov[1]=prot0.GetCov()[4];//sigmayz
1398 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
1399 if (!track0->Update(pointPos,pointCov)) break;
1401 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
1402 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
1403 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
1405 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
1406 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
1407 pointCov[0]=prot1.GetCov()[3];//simay^2
1408 pointCov[1]=prot1.GetCov()[4];//sigmayz
1409 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
1410 if (!track1->Update(pointPos,pointCov)) break;
1414 if (npoints2<npoints) return 0;
1415 AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
1416 track1->Rotate(track0->GetAlpha());
1417 AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
1419 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
1420 "point0.="<<&pointArray0<< // points
1421 "point1.="<<&pointArray1<< // distorted points
1422 "trackIn.="<<&track<< // original track
1423 "track0.="<<track0<< // fitted track
1424 "track1.="<<track1<< // fitted distorted track
1426 new(&trackIn) AliExternalTrackParam(*track0);
1435 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
1437 // create the distortion tree on a mesh with granularity given by step
1438 // return the tree with distortions at given position
1439 // Map is created on the mesh with given step size
1441 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
1443 for (Double_t x= -250; x<250; x+=step){
1444 for (Double_t y= -250; y<250; y+=step){
1445 Double_t r = TMath::Sqrt(x*x+y*y);
1447 if (r>250) continue;
1448 for (Double_t z= -250; z<250; z+=step){
1449 Int_t roc=(z>0)?0:18;
1453 Double_t phi = TMath::ATan2(y,x);
1454 DistortPoint(xyz,roc);
1455 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1456 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
1457 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
1458 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
1459 Double_t dx = xyz[0]-x;
1460 Double_t dy = xyz[1]-y;
1461 Double_t dz = xyz[2]-z;
1463 Double_t drphi=(phi1-phi)*r;
1464 (*pcstream)<<"distortion"<<
1465 "x="<<x<< // original position
1470 "x1="<<xyz[0]<< // distorted position
1476 "dx="<<dx<< // delta position
1486 TFile f(Form("correction%s.root",GetName()));
1487 TTree * tree = (TTree*)f.Get("distortion");
1488 TTree * tree2= tree->CopyTree("1");
1489 tree2->SetName(Form("dist%s",GetName()));
1490 tree2->SetDirectory(0);
1498 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
1501 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1502 // calculates partial distortions
1503 // Partial distortion is stored in the resulting tree
1504 // Output is storred in the file distortion_<dettype>_<partype>.root
1505 // Partial distortion is stored with the name given by correction name
1508 // Parameters of function:
1509 // input - input tree
1510 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex
1511 // ppype - parameter type
1512 // corrArray - array with partial corrections
1513 // step - skipe entries - if 1 all entries processed - it is slow
1514 // debug 0 if debug on also space points dumped - it is slow
1516 const Double_t kMaxSnp = 0.85;
1517 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1518 // const Double_t kB2C=-0.299792458e-3;
1519 const Int_t kMinEntries=50;
1520 Double_t phi,theta, snp, mean,rms, entries;
1521 tinput->SetBranchAddress("theta",&theta);
1522 tinput->SetBranchAddress("phi", &phi);
1523 tinput->SetBranchAddress("snp",&snp);
1524 tinput->SetBranchAddress("mean",&mean);
1525 tinput->SetBranchAddress("rms",&rms);
1526 tinput->SetBranchAddress("entries",&entries);
1527 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
1529 Int_t nentries=tinput->GetEntries();
1530 Int_t ncorr=corrArray->GetEntries();
1531 Double_t corrections[100]={0}; //
1533 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1536 if (dtype==0) {refX=85.; dir=-1;}
1537 if (dtype==1) {refX=275.; dir=1;}
1538 if (dtype==2) {refX=85.; dir=-1;}
1539 if (dtype==3) {refX=360.; dir=-1;}
1541 for (Int_t ientry=0; ientry<nentries; ientry+=step){
1542 tinput->GetEntry(ientry);
1543 if (TMath::Abs(snp)>kMaxSnp) continue;
1548 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
1549 Double_t bz=AliTrackerBase::GetBz();
1550 if (refX>10. && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
1551 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
1552 AliExternalTrackParam track(refX,phi,tPar,cov);
1556 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1557 if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature
1558 (*pcstream)<<"fit"<<
1559 "bz="<<bz<< // magnetic filed used
1560 "dtype="<<dtype<< // detector match type
1561 "ptype="<<ptype<< // parameter type
1562 "theta="<<theta<< // theta
1563 "phi="<<phi<< // phi
1564 "snp="<<snp<< // snp
1565 "mean="<<mean<< // mean dist value
1566 "rms="<<rms<< // rms
1567 "gx="<<xyz[0]<< // global position at reference
1568 "gy="<<xyz[1]<< // global position at reference
1569 "gz="<<xyz[2]<< // global position at reference
1570 "dRrec="<<dRrec<< // delta Radius in reconstruction
1571 "id="<<id<< // track id
1572 "entries="<<entries;// number of entries in bin
1574 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1575 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1576 corrections[icorr]=0;
1577 if (entries>kMinEntries){
1578 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1579 AliExternalTrackParam *trackOut = 0;
1580 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1581 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1582 if (dtype==0) {refX=85.; dir=-1;}
1583 if (dtype==1) {refX=275.; dir=1;}
1584 if (dtype==2) {refX=0; dir=-1;}
1585 if (dtype==3) {refX=360.; dir=-1;}
1588 AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1589 trackOut->Rotate(trackIn.GetAlpha());
1590 trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1592 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1595 corrections[icorr]=0;
1597 if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature
1600 (*pcstream)<<"fit"<<
1601 Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
1602 Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
1603 // for points it is neccessary
1605 (*pcstream)<<"fit"<<"\n";
1612 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1614 // Make a laser fit tree for global minimization
1616 const Double_t cutErrY=0.1;
1617 const Double_t cutErrZ=0.1;
1618 const Double_t kEpsilon=0.00000001;
1623 AliTPCLaserTrack *ltr=0;
1624 AliTPCLaserTrack::LoadTracks();
1625 tree->SetBranchAddress("dY.",&vecdY);
1626 tree->SetBranchAddress("dZ.",&vecdZ);
1627 tree->SetBranchAddress("eY.",&veceY);
1628 tree->SetBranchAddress("eZ.",&veceZ);
1629 tree->SetBranchAddress("LTr.",<r);
1630 Int_t entries= tree->GetEntries();
1631 TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1632 Double_t bz=AliTrackerBase::GetBz();
1635 for (Int_t ientry=0; ientry<entries; ientry++){
1636 tree->GetEntry(ientry);
1637 if (!ltr->GetVecGX()){
1638 ltr->UpdatePoints();
1640 TVectorD * delta= (itype==0)? vecdY:vecdZ;
1641 TVectorD * err= (itype==0)? veceY:veceZ;
1643 for (Int_t irow=0; irow<159; irow++){
1644 Int_t nentries = 1000;
1645 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1646 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1648 Double_t phi =(*ltr->GetVecPhi())[irow];
1649 Double_t theta =ltr->GetTgl();
1650 Double_t mean=delta->GetMatrixArray()[irow];
1651 Double_t gx=0,gy=0,gz=0;
1652 Double_t snp = (*ltr->GetVecP2())[irow];
1653 Double_t rms = 0.1+err->GetMatrixArray()[irow];
1654 gx = (*ltr->GetVecGX())[irow];
1655 gy = (*ltr->GetVecGY())[irow];
1656 gz = (*ltr->GetVecGZ())[irow];
1657 Int_t bundle= ltr->GetBundle();
1660 // get delta R used in reconstruction
1661 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
1662 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1663 const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1664 Double_t xyz0[3]={gx,gy,gz};
1665 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1667 // old ExB correction
1669 if(recoParam&&recoParam->GetUseExBCorrection()) {
1670 Double_t xyz1[3]={gx,gy,gz};
1671 calib->GetExB()->Correct(xyz0,xyz1);
1672 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1675 if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1676 Float_t xyz1[3]={gx,gy,gz};
1677 Int_t sector=(gz>0)?0:18;
1678 correction->CorrectPoint(xyz1, sector);
1679 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1684 (*pcstream)<<"fit"<<
1685 "bz="<<bz<< // magnetic filed used
1686 "dtype="<<dtype<< // detector match type
1687 "ptype="<<itype<< // parameter type
1688 "theta="<<theta<< // theta
1689 "phi="<<phi<< // phi
1690 "snp="<<snp<< // snp
1691 "mean="<<mean<< // mean dist value
1692 "rms="<<rms<< // rms
1693 "gx="<<gx<< // global position
1694 "gy="<<gy<< // global position
1695 "gz="<<gz<< // global position
1696 "dRrec="<<dRrec<< // delta Radius in reconstruction
1697 "id="<<bundle<< //bundle
1698 "entries="<<nentries;// number of entries in bin
1701 Double_t ky = TMath::Tan(TMath::ASin(snp));
1702 Int_t ncorr = corrArray->GetEntries();
1703 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
1704 Double_t phi0 = TMath::ATan2(gy,gx);
1705 Double_t distortions[1000]={0};
1706 Double_t distortionsR[1000]={0};
1707 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1708 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1709 Float_t distPoint[3]={gx,gy,gz};
1710 Int_t sector= (gz>0)? 0:18;
1712 corr->DistortPoint(distPoint, sector);
1714 // Double_t value=distPoint[2]-gz;
1716 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1717 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1718 Double_t drphi= r0*(phi1-phi0);
1719 Double_t dr = r1-r0;
1720 distortions[icorr] = drphi-ky*dr;
1721 distortionsR[icorr] = dr;
1723 (*pcstream)<<"fit"<<
1724 Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
1725 Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
1727 (*pcstream)<<"fit"<<"\n";
1735 void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
1737 // make a distortion map out ou fthe residual histogram
1738 // Results are written to the debug streamer - pcstream
1740 // his0 - input (4D) residual histogram
1741 // pcstream - file to write the tree
1743 // marian.ivanov@cern.ch
1744 const Int_t kMinEntries=50;
1745 Int_t nbins1=his0->GetAxis(1)->GetNbins();
1746 Int_t first1=his0->GetAxis(1)->GetFirst();
1747 Int_t last1 =his0->GetAxis(1)->GetLast();
1749 Double_t bz=AliTrackerBase::GetBz();
1750 Int_t idim[4]={0,1,2,3};
1751 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
1753 his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1754 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1755 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
1756 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
1757 Int_t first3 = his1->GetAxis(3)->GetFirst();
1758 Int_t last3 = his1->GetAxis(3)->GetLast();
1761 for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){ // axis 3 - local angle
1762 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1763 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1765 his1->GetAxis(3)->SetRangeUser(-1,1);
1768 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
1769 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1770 Int_t first2 = his3->GetAxis(2)->GetFirst();
1771 Int_t last2 = his3->GetAxis(2)->GetLast();
1773 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1774 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1775 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1776 TH1 * hisDelta = his3->Projection(0);
1778 Double_t entries = hisDelta->GetEntries();
1779 Double_t mean=0, rms=0;
1780 if (entries>kMinEntries){
1781 mean = hisDelta->GetMean();
1782 rms = hisDelta->GetRMS();
1784 (*pcstream)<<hname<<
1790 "entries="<<entries<<
1795 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1807 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1809 // Store object in the OCDB
1810 // By default the object is stored in the current directory
1811 // default comment consit of user name and the date
1813 TString ocdbStorage="";
1814 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1815 AliCDBMetaData *metaData= new AliCDBMetaData();
1816 metaData->SetObjectClassName("AliTPCCorrection");
1817 metaData->SetResponsible("Marian Ivanov");
1818 metaData->SetBeamPeriod(1);
1819 metaData->SetAliRootVersion("05-25-01"); //root version
1820 TString userName=gSystem->GetFromPipe("echo $USER");
1821 TString date=gSystem->GetFromPipe("date");
1823 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1824 if (comment) metaData->SetComment(comment);
1826 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1827 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1828 gStorage->Put(this, (*id1), metaData);
1832 void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
1834 // Fast method to simulate the influence of the given distortion on the vertex reconstruction
1837 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
1838 if (!magF) AliError("Magneticd field - not initialized");
1839 Double_t bz = magF->SolenoidField(); //field in kGauss
1840 printf("bz: %lf\n",bz);
1841 AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
1843 TObjArray aTrk; // Original Track array of Aside
1844 TObjArray daTrk; // Distorted Track array of A side
1845 UShort_t *aId = new UShort_t[nTracks]; // A side Track ID
1848 UShort_t *cId = new UShort_t [nTracks];
1850 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1851 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
1852 fpt.SetParameters(7.24,0.120);
1854 for(Int_t nt=0; nt<nTracks; nt++){
1855 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
1856 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
1857 Double_t pt = fpt.GetRandom(); // momentum for f1
1858 // printf("phi %lf eta %lf pt %lf\n",phi,eta,pt);
1860 if(gRandom->Rndm() < 0.5){
1866 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1868 pxyz[0]=pt*TMath::Cos(phi);
1869 pxyz[1]=pt*TMath::Sin(phi);
1870 pxyz[2]=pt*TMath::Tan(theta);
1871 Double_t cv[21]={0};
1872 AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
1876 AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir, NULL);
1878 if (pcstream) (*pcstream)<<"track"<<
1884 if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
1888 Int_t nn=aTrk.GetEntriesFast();
1891 }else if(( eta<-0.07 )&&( eta>-etaCuts )){
1895 Int_t nn=cTrk.GetEntriesFast();
1900 }// end of track loop
1902 vertexer->SetTPCMode();
1903 vertexer->SetConstraintOff();
1905 aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));
1906 avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
1907 cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));
1908 cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
1909 if (pcstream) (*pcstream)<<"vertex"<<
1910 "x="<<orgVertex[0]<<
1911 "y="<<orgVertex[1]<<
1912 "z="<<orgVertex[2]<<
1913 "av.="<<&aV<< // distorted vertex A side
1914 "cv.="<<&cV<< // distroted vertex C side
1915 "avO.="<<&avOrg<< // original vertex A side
1922 void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
1924 // make correction available for visualization using
1925 // TFormula, TFX and TTree::Draw
1926 // important in order to check corrections and also compute dervied variables
1927 // e.g correction partial derivatives
1929 // NOTE - class is not owner of correction
1931 if (!fgVisualCorrection) fgVisualCorrection=new TObjArray;
1932 if (position>=fgVisualCorrection->GetEntriesFast()) fgVisualCorrection->Expand(position*2);
1933 fgVisualCorrection->AddAt(corr, position);
1938 Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t r, Double_t kZ, Int_t axisType, Int_t corrType){
1940 // calculate the correction at given position - check the geffCorr
1942 if (!fgVisualCorrection) return 0;
1943 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1944 if (!corr) return 0;
1946 Double_t phi=sector*TMath::Pi()/9.;
1947 Double_t gx = r*TMath::Cos(phi);
1948 Double_t gy = r*TMath::Sin(phi);
1950 Int_t nsector=(gz>0) ? 0:18;
1954 Float_t distPoint[3]={gx,gy,gz};
1955 corr->DistortPoint(distPoint, nsector);
1956 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1957 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1958 Double_t phi0=TMath::ATan2(gy,gx);
1959 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1960 if (axisType==0) return r1-r0;
1961 if (axisType==1) return (phi1-phi0)*r0;
1962 if (axisType==2) return distPoint[2]-gz;
1966 Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
1968 // return correction at given x,y,z
1970 if (!fgVisualCorrection) return 0;
1971 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1972 if (!corr) return 0;
1973 Double_t phi0= TMath::ATan2(gy,gx);
1974 Int_t nsector=(gz>0) ? 0:18;
1975 Float_t distPoint[3]={gx,gy,gz};
1976 corr->DistortPoint(distPoint, nsector);
1977 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1978 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1979 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1980 if (axisType==0) return r1-r0;
1981 if (axisType==1) return (phi1-phi0)*r0;
1982 if (axisType==2) return distPoint[2]-gz;