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.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
92 const Double_t AliTPCCorrection::fgkOFCRadius= 254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
93 const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
94 const Double_t AliTPCCorrection::fgkCathodeV = -100000.0; // Cathode Voltage (volts)
95 const Double_t AliTPCCorrection::fgkGG = -70.0; // Gating Grid voltage (volts)
97 const Double_t AliTPCCorrection::fgkdvdE = 0.0024; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
99 const Double_t AliTPCCorrection::fgkEM = -1.602176487e-19/9.10938215e-31; // charge/mass in [C/kg]
100 const Double_t AliTPCCorrection::fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)]
102 // FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
103 const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
104 83.06, 83.5, 84.0, 84.5, 85.0, 85.5, 86.0, 86.5,
105 87.0, 87.5, 88.0, 88.5, 89.0, 89.5, 90.0, 90.5, 91.0, 92.0,
106 93.0, 94.0, 95.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
107 110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
108 130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
109 150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
110 170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
111 190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
112 210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
113 230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 241.0, 242.0, 243.0, 244.0,
114 245.0, 245.5, 246.0, 246.5, 247.0, 247.5, 248.0, 248.5, 249.0, 249.5,
115 250.0, 250.5, 251.0, 251.5, 252.0, 252.5, 253.0, 253.5, 254.0, 254.5
118 const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = {
119 -249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
120 -240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
121 -220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
122 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
123 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
124 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
125 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
126 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
127 -100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0,
128 -80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0,
129 -60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0,
130 -40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0,
131 -20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0,
132 -1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0,
133 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
134 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0,
135 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0,
136 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0,
137 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0,
138 102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
139 122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
140 142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
141 162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
142 182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
143 202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0,
144 222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0,
145 242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ;
149 AliTPCCorrection::AliTPCCorrection()
150 : TNamed("correction_unity","unity"),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
153 // default constructor
155 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
157 // Initialization of interpolation points
158 for (Int_t i = 0; i<kNPhi; i++) {
159 fgkPhiList[i] = TMath::TwoPi()*i/(kNPhi-1);
164 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
165 : TNamed(name,title),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
168 // default constructor, that set the name and title
170 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
172 // Initialization of interpolation points
173 for (Int_t i = 0; i<kNPhi; i++) {
174 fgkPhiList[i] = TMath::TwoPi()*i/(kNPhi-1);
179 AliTPCCorrection::~AliTPCCorrection() {
181 // virtual destructor
185 void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
187 // Corrects the initial coordinates x (cartesian coordinates)
188 // according to the given effect (inherited classes)
189 // roc represents the TPC read out chamber (offline numbering convention)
192 GetCorrection(x,roc,dx);
193 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
196 void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
198 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
199 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
200 // roc represents the TPC read out chamber (offline numbering convention)
203 GetCorrection(x,roc,dx);
204 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
207 void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
209 // Distorts the initial coordinates x (cartesian coordinates)
210 // according to the given effect (inherited classes)
211 // roc represents the TPC read out chamber (offline numbering convention)
214 GetDistortion(x,roc,dx);
215 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
218 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
220 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
221 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
222 // roc represents the TPC read out chamber (offline numbering convention)
225 GetDistortion(x,roc,dx);
226 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
229 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
231 // This function delivers the correction values dx in respect to the inital coordinates x
232 // roc represents the TPC read out chamber (offline numbering convention)
233 // Note: The dx is overwritten by the inherited effectice class ...
235 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
238 void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
240 // This function delivers the distortion values dx in respect to the inital coordinates x
241 // roc represents the TPC read out chamber (offline numbering convention)
243 GetCorrection(x,roc,dx);
244 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
247 void AliTPCCorrection::Init() {
249 // Initialization funtion (not used at the moment)
253 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
259 void AliTPCCorrection::Print(Option_t* /*option*/) const {
261 // Print function to check which correction classes are used
262 // option=="d" prints details regarding the setted magnitude
263 // option=="a" prints the C0 and C1 coefficents for calibration purposes
265 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
268 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
270 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
271 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
276 //SetOmegaTauT1T2(omegaTau, t1, t2);
279 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
281 // Simple plot functionality.
282 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
283 // in respect to position z within the XY plane.
284 // The histogramm has nx times ny entries.
286 AliTPCParam* tpcparam = new AliTPCParamSR;
288 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
289 nx,-250.,250.,ny,-250.,250.);
292 Int_t roc=z>0.?0:18; // FIXME
293 for (Int_t iy=1;iy<=ny;++iy) {
294 x[1]=h->GetYaxis()->GetBinCenter(iy);
295 for (Int_t ix=1;ix<=nx;++ix) {
296 x[0]=h->GetXaxis()->GetBinCenter(ix);
297 GetCorrection(x,roc,dx);
298 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
299 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
300 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
301 h->SetBinContent(ix,iy,r1-r0);
304 h->SetBinContent(ix,iy,0.);
311 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
313 // Simple plot functionality.
314 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
315 // in respect to position z within the XY plane.
316 // The histogramm has nx times ny entries.
319 AliTPCParam* tpcparam = new AliTPCParamSR;
321 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
322 nx,-250.,250.,ny,-250.,250.);
325 Int_t roc=z>0.?0:18; // FIXME
326 for (Int_t iy=1;iy<=ny;++iy) {
327 x[1]=h->GetYaxis()->GetBinCenter(iy);
328 for (Int_t ix=1;ix<=nx;++ix) {
329 x[0]=h->GetXaxis()->GetBinCenter(ix);
330 GetCorrection(x,roc,dx);
331 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
332 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
333 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
334 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
336 Float_t dphi=phi1-phi0;
337 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
338 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
340 h->SetBinContent(ix,iy,r0*dphi);
343 h->SetBinContent(ix,iy,0.);
350 TH2F* AliTPCCorrection::CreateHistoDZinXY(Float_t z,Int_t nx,Int_t ny) {
352 // Simple plot functionality.
353 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
354 // in respect to position z within the XY plane.
355 // The histogramm has nx times ny entries.
358 AliTPCParam* tpcparam = new AliTPCParamSR;
360 TH2F *h=CreateTH2F("dz_xy",GetTitle(),"x [cm]","y [cm]","dz [cm]",
361 nx,-250.,250.,ny,-250.,250.);
364 Int_t roc=z>0.?0:18; // FIXME
365 for (Int_t iy=1;iy<=ny;++iy) {
366 x[1]=h->GetYaxis()->GetBinCenter(iy);
367 for (Int_t ix=1;ix<=nx;++ix) {
368 x[0]=h->GetXaxis()->GetBinCenter(ix);
369 GetCorrection(x,roc,dx);
370 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
371 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
372 h->SetBinContent(ix,iy,dx[2]);
375 h->SetBinContent(ix,iy,0.);
382 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
384 // Simple plot functionality.
385 // Returns a 2d hisogram which represents the corrections in r direction (dr)
386 // in respect to angle phi within the ZR plane.
387 // The histogramm has nx times ny entries.
389 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
390 nz,-250.,250.,nr,85.,250.);
392 for (Int_t ir=1;ir<=nr;++ir) {
393 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
394 x[0]=radius*TMath::Cos(phi);
395 x[1]=radius*TMath::Sin(phi);
396 for (Int_t iz=1;iz<=nz;++iz) {
397 x[2]=h->GetXaxis()->GetBinCenter(iz);
398 Int_t roc=x[2]>0.?0:18; // FIXME
399 GetCorrection(x,roc,dx);
400 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
401 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
402 h->SetBinContent(iz,ir,r1-r0);
409 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
411 // Simple plot functionality.
412 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
413 // in respect to angle phi within the ZR plane.
414 // The histogramm has nx times ny entries.
416 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
417 nz,-250.,250.,nr,85.,250.);
419 for (Int_t iz=1;iz<=nz;++iz) {
420 x[2]=h->GetXaxis()->GetBinCenter(iz);
421 Int_t roc=x[2]>0.?0:18; // FIXME
422 for (Int_t ir=1;ir<=nr;++ir) {
423 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
424 x[0]=radius*TMath::Cos(phi);
425 x[1]=radius*TMath::Sin(phi);
426 GetCorrection(x,roc,dx);
427 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
428 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
429 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
431 Float_t dphi=phi1-phi0;
432 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
433 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
435 h->SetBinContent(iz,ir,r0*dphi);
441 TH2F* AliTPCCorrection::CreateHistoDZinZR(Float_t phi,Int_t nz,Int_t nr) {
443 // Simple plot functionality.
444 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
445 // in respect to angle phi within the ZR plane.
446 // The histogramm has nx times ny entries.
448 TH2F *h=CreateTH2F("dz_zr",GetTitle(),"z [cm]","r [cm]","dz [cm]",
449 nz,-250.,250.,nr,85.,250.);
451 for (Int_t ir=1;ir<=nr;++ir) {
452 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
453 x[0]=radius*TMath::Cos(phi);
454 x[1]=radius*TMath::Sin(phi);
455 for (Int_t iz=1;iz<=nz;++iz) {
456 x[2]=h->GetXaxis()->GetBinCenter(iz);
457 Int_t roc=x[2]>0.?0:18; // FIXME
458 GetCorrection(x,roc,dx);
459 h->SetBinContent(iz,ir,dx[2]);
467 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
468 const char *xlabel,const char *ylabel,const char *zlabel,
469 Int_t nbinsx,Double_t xlow,Double_t xup,
470 Int_t nbinsy,Double_t ylow,Double_t yup) {
472 // Helper function to create a 2d histogramm of given size
478 while (gDirectory->FindObject(hname.Data())) {
485 TH2F *h=new TH2F(hname.Data(),title,
488 h->GetXaxis()->SetTitle(xlabel);
489 h->GetYaxis()->SetTitle(ylabel);
490 h->GetZaxis()->SetTitle(zlabel);
495 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
497 void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
498 const Double_t er[kNZ][kNR], Double_t &erValue ) {
500 // Interpolate table - 2D interpolation
502 Double_t saveEr[10] ;
504 Search( kNZ, fgkZList, z, fJLow ) ;
505 Search( kNR, fgkRList, r, fKLow ) ;
506 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
507 if ( fKLow < 0 ) fKLow = 0 ;
508 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
509 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
511 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
512 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
514 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
518 void AliTPCCorrection::Interpolate3DEdistortion( const Int_t order, const Double_t r, const Float_t phi, const Double_t z,
519 const Double_t er[kNZ][kNPhi][kNR], const Double_t ephi[kNZ][kNPhi][kNR], const Double_t ez[kNZ][kNPhi][kNR],
520 Double_t &erValue, Double_t &ephiValue, Double_t &ezValue) {
522 // Interpolate table - 3D interpolation
525 Double_t saveEr[10], savedEr[10] ;
526 Double_t saveEphi[10], savedEphi[10] ;
527 Double_t saveEz[10], savedEz[10] ;
529 Search( kNZ, fgkZList, z, fILow ) ;
530 Search( kNPhi, fgkPhiList, z, fJLow ) ;
531 Search( kNR, fgkRList, r, fKLow ) ;
533 if ( fILow < 0 ) fILow = 0 ; // check if out of range
534 if ( fJLow < 0 ) fJLow = 0 ;
535 if ( fKLow < 0 ) fKLow = 0 ;
537 if ( fILow + order >= kNZ - 1 ) fILow = kNZ - 1 - order ;
538 if ( fJLow + order >= kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
539 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
541 for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
542 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
543 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r ) ;
544 saveEphi[j-fJLow] = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
545 saveEz[j-fJLow] = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r ) ;
547 savedEr[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi ) ;
548 savedEphi[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
549 savedEz[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi ) ;
551 erValue = Interpolate( &fgkZList[fILow], savedEr, order, z ) ;
552 ephiValue = Interpolate( &fgkZList[fILow], savedEphi, order, z ) ;
553 ezValue = Interpolate( &fgkZList[fILow], savedEz, order, z ) ;
557 Double_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
558 const Int_t nx, const Int_t ny, const Double_t xv[], const Double_t yv[],
559 const TMatrixD &array ) {
561 // Interpolate table (TMatrix format) - 2D interpolation
564 static Int_t jlow = 0, klow = 0 ;
565 Double_t saveArray[10] ;
567 Search( nx, xv, x, jlow ) ;
568 Search( ny, yv, y, klow ) ;
569 if ( jlow < 0 ) jlow = 0 ; // check if out of range
570 if ( klow < 0 ) klow = 0 ;
571 if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
572 if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
574 for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
576 Double_t *ajkl = &((TMatrixD&)array)(j,klow);
577 saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
580 return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
584 Double_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x, const Double_t y, const Double_t z,
585 const Int_t nx, const Int_t ny, const Int_t nz,
586 const Double_t xv[], const Double_t yv[], const Double_t zv[],
587 TMatrixD **arrayofArrays ) {
589 // Interpolate table (TMatrix format) - 3D interpolation
592 static Int_t ilow = 0, jlow = 0, klow = 0 ;
593 Double_t saveArray[10], savedArray[10] ;
595 Search( nx, xv, x, ilow ) ;
596 Search( ny, yv, y, jlow ) ;
597 Search( nz, zv, z, klow ) ;
599 if ( ilow < 0 ) ilow = 0 ; // check if out of range
600 if ( jlow < 0 ) jlow = 0 ;
601 if ( klow < 0 ) klow = 0 ;
603 if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
604 if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
605 if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
607 for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
609 TMatrixD &table = *arrayofArrays[k] ;
610 for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
612 saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
614 savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
616 return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
622 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
623 const Int_t order, const Double_t x ) {
625 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
629 if ( order == 2 ) { // Quadratic Interpolation = 2
630 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
631 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
632 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
633 } else { // Linear Interpolation = 1
634 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
642 void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
644 // Search an ordered table by starting at the most recently used point
647 Long_t middle, high ;
648 Int_t ascend = 0, increment = 1 ;
650 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
652 if ( low < 0 || low > n-1 ) {
653 low = -1 ; high = n ;
654 } else { // Ordered Search phase
655 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
656 if ( low == n-1 ) return ;
658 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
661 high = low + increment ;
662 if ( high > n-1 ) { high = n ; break ; }
665 if ( low == 0 ) { low = -1 ; return ; }
667 while ( (Int_t)( x < xArray[low] ) == ascend ) {
670 if ( increment >= high ) { low = -1 ; break ; }
671 else low = high - increment ;
676 while ( (high-low) != 1 ) { // Binary Search Phase
677 middle = ( high + low ) / 2 ;
678 if ( (Int_t)( x >= xArray[middle] ) == ascend )
684 if ( x == xArray[n-1] ) low = n-2 ;
685 if ( x == xArray[0] ) low = 0 ;
689 void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
690 TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
691 const Int_t rows, const Int_t columns, const Int_t iterations,
692 const Bool_t rocDisplacement ) {
694 // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
696 // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
697 // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
698 // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
699 // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
700 // you wish to solve Laplaces equation however it should not contain random numbers or you will get
701 // random numbers back as a solution.
702 // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
703 // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
704 // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
705 // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
706 // number of points and solves the problem again. This happens several times until the maximum number
707 // of points has been included in the array.
709 // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
710 // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
712 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
714 // Original code by Jim Thomas (STAR TPC Collaboration)
717 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
719 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
720 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
721 const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
723 TMatrixD arrayEr(rows,columns) ;
724 TMatrixD arrayEz(rows,columns) ;
726 //Check that number of rows and columns is suitable for a binary expansion
728 if ( !IsPowerOfTwo(rows-1) ) {
729 AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
732 if ( !IsPowerOfTwo(columns-1) ) {
733 AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
737 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
738 // Allow for different size grid spacing in R and Z directions
739 // Use a binary expansion of the size of the matrix to speed up the solution of the problem
741 Int_t iOne = (rows-1)/4 ;
742 Int_t jOne = (columns-1)/4 ;
743 // Solve for N in 2**N, add one.
744 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
746 for ( Int_t count = 0 ; count < loops ; count++ ) {
747 // Loop while the matrix expands & the resolution increases.
749 Float_t tempGridSizeR = gridSizeR * iOne ;
750 Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
751 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
753 // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
754 std::vector<float> coef1(rows) ;
755 std::vector<float> coef2(rows) ;
757 for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
758 Float_t radius = fgkIFCRadius + i*gridSizeR ;
759 coef1[i] = 1.0 + tempGridSizeR/(2*radius);
760 coef2[i] = 1.0 - tempGridSizeR/(2*radius);
763 TMatrixD sumChargeDensity(rows,columns) ;
765 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
766 Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
767 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
768 if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
770 // Add up all enclosed charge density contributions within 1/2 unit in all directions
771 Float_t weight = 0.0 ;
773 sumChargeDensity(i,j) = 0.0 ;
774 for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
775 for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
776 if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
779 // Note that this is cylindrical geometry
780 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
781 sum += weight*radius ;
784 sumChargeDensity(i,j) /= sum ;
786 sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
790 for ( Int_t k = 1 ; k <= iterations; k++ ) {
791 // Solve Poisson's Equation
792 // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
793 // as interations increase.
794 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
795 Float_t overRelaxM1 = overRelax - 1.0 ;
796 Float_t overRelaxtempFourth, overRelaxcoef5 ;
797 overRelaxtempFourth = overRelax * tempFourth ;
798 overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
800 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
801 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
803 arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
804 + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
805 - overRelaxcoef5 * arrayV(i,j)
806 + coef1[i] * arrayV(i+iOne,j)
807 + sumChargeDensity(i,j)
808 ) * overRelaxtempFourth;
812 if ( k == iterations ) {
813 // After full solution is achieved, copy low resolution solution into higher res array
814 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
815 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
818 arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
819 if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
822 arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
823 if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
825 if ( iOne > 1 && jOne > 1 ) {
826 arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
827 if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
828 if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
829 // Note that this leaves a point at the upper left and lower right corners uninitialized.
830 // -> Not a big deal.
839 iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
840 jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
842 sumChargeDensity.Clear();
845 // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
846 for ( Int_t j = 0 ; j < columns ; j++ ) {
847 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
848 arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
849 arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
852 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
853 for ( Int_t i = 0 ; i < rows ; i++) {
854 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
855 arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
856 arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
859 for ( Int_t i = 0 ; i < rows ; i++) {
860 // Note: go back and compare to old version of this code. See notes below.
861 // JT Test ... attempt to divide by real Ez not Ez to first order
862 for ( Int_t j = 0 ; j < columns ; j++ ) {
863 arrayEz(i,j) += ezField;
864 // This adds back the overall Z gradient of the field (main E field component)
866 // Warning: (-=) assumes you are using an error potetial without the overall Field included
869 // Integrate Er/Ez from Z to zero
870 for ( Int_t j = 0 ; j < columns ; j++ ) {
871 for ( Int_t i = 0 ; i < rows ; i++ ) {
873 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
874 arrayErOverEz(i,j) = 0.0 ;
875 arrayDeltaEz(i,j) = 0.0 ;
877 for ( Int_t k = j ; k < columns ; k++ ) {
878 arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
879 arrayDeltaEz(i,j) += index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
880 if ( index != 4 ) index = 4; else index = 2 ;
883 arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
884 arrayDeltaEz(i,j) -= (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
887 arrayErOverEz(i,j) += (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
888 -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
889 arrayDeltaEz(i,j) += (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
890 -2.5*(arrayEz(i,columns-1)-ezField));
892 if ( j == columns-2 ) {
893 arrayErOverEz(i,j) = (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
894 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
895 arrayDeltaEz(i,j) = (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
896 +1.5*(arrayEz(i,columns-1)-ezField) ) ;
898 if ( j == columns-1 ) {
899 arrayErOverEz(i,j) = 0.0 ;
900 arrayDeltaEz(i,j) = 0.0 ;
905 // calculate z distortion from the integrated Delta Ez residuals
906 // and include the aquivalence (Volt to cm) of the ROC shift !!
908 for ( Int_t j = 0 ; j < columns ; j++ ) {
909 for ( Int_t i = 0 ; i < rows ; i++ ) {
911 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
912 arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
914 // ROC Potential in cm aquivalent
915 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
916 if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift; // add the ROC misaligment
926 void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
927 TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
928 const Int_t rows, const Int_t columns, const Int_t phislices,
929 const Float_t deltaphi, const Int_t iterations, const Int_t symmetry,
930 Bool_t rocDisplacement ) {
932 // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
934 // NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
935 // The number of rows and COLUMNS can be different.
938 // COLUMNS == 2**N + 1
939 // PHISLICES == Arbitrary but greater than 3
941 // DeltaPhi in Radians
943 // SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
944 // = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
946 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
948 const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
950 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
951 const Float_t gridSizePhi = deltaphi ;
952 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
953 const Float_t ratioPhi = gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
954 const Float_t ratioZ = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
956 TMatrixD arrayE(rows,columns) ;
958 // Check that the number of rows and columns is suitable for a binary expansion
959 if ( !IsPowerOfTwo((rows-1)) ) {
960 AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
962 if ( !IsPowerOfTwo((columns-1)) ) {
963 AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
965 if ( phislices <= 3 ) {
966 AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
968 if ( phislices > 1000 ) {
969 AliError("Poisson3D phislices > 1000 is not allowed (nor wise) ");
972 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
973 // Allow for different size grid spacing in R and Z directions
974 // Use a binary expansion of the matrix to speed up the solution of the problem
976 Int_t loops, mplus, mminus, signplus, signminus ;
977 Int_t ione = (rows-1)/4 ;
978 Int_t jone = (columns-1)/4 ;
979 loops = TMath::Max(ione, jone) ; // Calculate the number of loops for the binary expansion
980 loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ; // Solve for N in 2**N
982 TMatrixD* arrayofSumChargeDensities[1000] ; // Create temporary arrays to store low resolution charge arrays
984 for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
986 for ( Int_t count = 0 ; count < loops ; count++ ) { // START the master loop and do the binary expansion
988 Float_t tempgridSizeR = gridSizeR * ione ;
989 Float_t tempratioPhi = ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
990 Float_t tempratioZ = ratioZ * ione * ione / ( jone * jone ) ;
992 std::vector<float> coef1(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
993 std::vector<float> coef2(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
994 std::vector<float> coef3(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
995 std::vector<float> coef4(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
997 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
998 Float_t radius = fgkIFCRadius + i*gridSizeR ;
999 coef1[i] = 1.0 + tempgridSizeR/(2*radius);
1000 coef2[i] = 1.0 - tempgridSizeR/(2*radius);
1001 coef3[i] = tempratioPhi/(radius*radius);
1002 coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
1005 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1006 TMatrixD &chargeDensity = *arrayofChargeDensities[m] ;
1007 TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
1008 for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
1009 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1010 for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
1011 if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1012 else { // Add up all enclosed charge density contributions within 1/2 unit in all directions
1013 Float_t weight = 0.0 ;
1015 sumChargeDensity(i,j) = 0.0 ;
1016 for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
1017 for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
1018 if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
1021 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1022 sum += weight*radius ;
1025 sumChargeDensity(i,j) /= sum ;
1027 sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
1032 for ( Int_t k = 1 ; k <= iterations; k++ ) {
1034 // over-relaxation index, >= 1 but < 2
1035 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1036 Float_t overRelaxM1 = overRelax - 1.0 ;
1038 std::vector<float> overRelaxcoef4(rows) ; // Do this the standard C++ way to avoid gcc extensions
1039 std::vector<float> overRelaxcoef5(rows) ; // Do this the standard C++ way to avoid gcc extensions
1041 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1042 overRelaxcoef4[i] = overRelax * coef4[i] ;
1043 overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
1046 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1048 mplus = m + 1; signplus = 1 ;
1049 mminus = m - 1 ; signminus = 1 ;
1050 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1051 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1052 if ( mminus < 0 ) mminus = 1 ;
1054 else if (symmetry==-1) { // Anti-symmetry in phi
1055 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1056 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1058 else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1059 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1060 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1062 TMatrixD& arrayV = *arrayofArrayV[m] ;
1063 TMatrixD& arrayVP = *arrayofArrayV[mplus] ;
1064 TMatrixD& arrayVM = *arrayofArrayV[mminus] ;
1065 TMatrixD& sumChargeDensity = *arrayofSumChargeDensities[m] ;
1067 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1068 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1070 arrayV(i,j) = ( coef2[i] * arrayV(i-ione,j)
1071 + tempratioZ * ( arrayV(i,j-jone) + arrayV(i,j+jone) )
1072 - overRelaxcoef5[i] * arrayV(i,j)
1073 + coef1[i] * arrayV(i+ione,j)
1074 + coef3[i] * ( signplus*arrayVP(i,j) + signminus*arrayVM(i,j) )
1075 + sumChargeDensity(i,j)
1076 ) * overRelaxcoef4[i] ;
1077 // Note: over-relax the solution at each step. This speeds up the convergance.
1082 if ( k == iterations ) { // After full solution is achieved, copy low resolution solution into higher res array
1083 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1084 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1087 arrayV(i+ione/2,j) = ( arrayV(i+ione,j) + arrayV(i,j) ) / 2 ;
1088 if ( i == ione ) arrayV(i-ione/2,j) = ( arrayV(0,j) + arrayV(ione,j) ) / 2 ;
1091 arrayV(i,j+jone/2) = ( arrayV(i,j+jone) + arrayV(i,j) ) / 2 ;
1092 if ( j == jone ) arrayV(i,j-jone/2) = ( arrayV(i,0) + arrayV(i,jone) ) / 2 ;
1094 if ( ione > 1 && jone > 1 ) {
1095 arrayV(i+ione/2,j+jone/2) = ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
1096 if ( i == ione ) arrayV(i-ione/2,j-jone/2) = ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
1097 if ( j == jone ) arrayV(i-ione/2,j-jone/2) = ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
1098 // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
1107 ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
1108 jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
1112 //Differentiate V(r) and solve for E(r) using special equations for the first and last row
1113 //Integrate E(r)/E(z) from point of origin to pad plane
1115 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1116 TMatrixD& arrayV = *arrayofArrayV[m] ;
1117 TMatrixD& eroverEz = *arrayofEroverEz[m] ;
1119 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1121 // Differentiate in R
1122 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1123 arrayE(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1124 arrayE(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1126 for ( Int_t i = 0 ; i < rows ; i++ ) {
1127 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1128 eroverEz(i,j) = 0.0 ;
1129 for ( Int_t k = j ; k < columns ; k++ ) {
1131 eroverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1132 if ( index != 4 ) index = 4; else index = 2 ;
1134 if ( index == 4 ) eroverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1135 if ( index == 2 ) eroverEz(i,j) +=
1136 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1137 if ( j == columns-2 ) eroverEz(i,j) =
1138 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1139 if ( j == columns-1 ) eroverEz(i,j) = 0.0 ;
1142 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1143 // eroverEz.Draw("surf") ; } // JT test
1146 //Differentiate V(r) and solve for E(phi)
1147 //Integrate E(phi)/E(z) from point of origin to pad plane
1149 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1151 mplus = m + 1; signplus = 1 ;
1152 mminus = m - 1 ; signminus = 1 ;
1153 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1154 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1155 if ( mminus < 0 ) mminus = 1 ;
1157 else if (symmetry==-1) { // Anti-symmetry in phi
1158 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1159 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1161 else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
1162 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1163 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1165 TMatrixD &arrayVP = *arrayofArrayV[mplus] ;
1166 TMatrixD &arrayVM = *arrayofArrayV[mminus] ;
1167 TMatrixD &ePhioverEz = *arrayofEPhioverEz[m] ;
1168 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1169 // Differentiate in Phi
1170 for ( Int_t i = 0 ; i < rows ; i++ ) {
1171 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1172 arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
1175 for ( Int_t i = 0 ; i < rows ; i++ ) {
1176 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1177 ePhioverEz(i,j) = 0.0 ;
1178 for ( Int_t k = j ; k < columns ; k++ ) {
1180 ePhioverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1181 if ( index != 4 ) index = 4; else index = 2 ;
1183 if ( index == 4 ) ePhioverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1184 if ( index == 2 ) ePhioverEz(i,j) +=
1185 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1186 if ( j == columns-2 ) ePhioverEz(i,j) =
1187 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1188 if ( j == columns-1 ) ePhioverEz(i,j) = 0.0 ;
1191 // if ( m == 5 ) { TCanvas* c2 = new TCanvas("arrayE","arrayE",50,50,840,600) ; c2 -> cd() ;
1192 // arrayE.Draw("surf") ; } // JT test
1196 // Differentiate V(r) and solve for E(z) using special equations for the first and last row
1197 // Integrate (E(z)-Ezstd) from point of origin to pad plane
1199 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1200 TMatrixD& arrayV = *arrayofArrayV[m] ;
1201 TMatrixD& deltaEz = *arrayofDeltaEz[m] ;
1203 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1204 for ( Int_t i = 0 ; i < rows ; i++) {
1205 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1206 arrayE(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1207 arrayE(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1210 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1212 for ( Int_t i = 0 ; i < rows ; i++ ) {
1213 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1214 deltaEz(i,j) = 0.0 ;
1215 for ( Int_t k = j ; k < columns ; k++ ) {
1216 deltaEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k) ;
1217 if ( index != 4 ) index = 4; else index = 2 ;
1219 if ( index == 4 ) deltaEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1) ;
1220 if ( index == 2 ) deltaEz(i,j) +=
1221 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
1222 if ( j == columns-2 ) deltaEz(i,j) =
1223 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
1224 if ( j == columns-1 ) deltaEz(i,j) = 0.0 ;
1227 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1228 // eroverEz.Draw("surf") ; } // JT test
1230 // calculate z distortion from the integrated Delta Ez residuals
1231 // and include the aquivalence (Volt to cm) of the ROC shift !!
1233 for ( Int_t j = 0 ; j < columns ; j++ ) {
1234 for ( Int_t i = 0 ; i < rows ; i++ ) {
1236 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1237 deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
1239 // ROC Potential in cm aquivalent
1240 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
1241 if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift; // add the ROC misaligment
1246 } // end loop over phi
1250 for ( Int_t k = 0 ; k < phislices ; k++ )
1252 arrayofSumChargeDensities[k]->Delete() ;
1261 Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1263 // Helperfunction: Check if integer is a power of 2
1266 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
1267 if ( j == 1 ) return(1) ; // True
1268 return(0) ; // False
1272 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
1274 // Fit the track parameters - without and with distortion
1275 // 1. Space points in the TPC are simulated along the trajectory
1276 // 2. Space points distorted
1277 // 3. Fits the non distorted and distroted track to the reference plane at refX
1278 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
1280 // trackIn - input track parameters
1281 // refX - reference X to fit the track
1282 // dir - direction - out=1 or in=-1
1283 // pcstream - debug streamer to check the results
1285 // see AliExternalTrackParam.h documentation:
1286 // track1.fP[0] - local y (rphi)
1288 // track1.fP[2] - sinus of local inclination angle
1289 // track1.fP[3] - tangent of deep angle
1290 // track1.fP[4] - 1/pt
1292 AliTPCROC * roc = AliTPCROC::Instance();
1293 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
1294 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
1295 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1296 const Double_t kMaxSnp = 0.85;
1297 const Double_t kSigmaY=0.1;
1298 const Double_t kSigmaZ=0.1;
1299 const Double_t kMaxR=500;
1300 const Double_t kMaxZ=500;
1301 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1305 AliExternalTrackParam track(trackIn); //
1307 AliTrackPointArray pointArray0(npoints0);
1308 AliTrackPointArray pointArray1(npoints0);
1310 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
1312 // simulate the track
1314 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
1315 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
1316 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
1318 xyz[0]+=gRandom->Gaus(0,0.00005);
1319 xyz[1]+=gRandom->Gaus(0,0.00005);
1320 xyz[2]+=gRandom->Gaus(0,0.00005);
1321 if (TMath::Abs(track.GetZ())>kMaxZ) break;
1322 if (TMath::Abs(track.GetX())>kMaxR) break;
1323 AliTrackPoint pIn0; // space point
1325 Int_t sector= (xyz[2]>0)? 0:18;
1326 pointArray0.GetPoint(pIn0,npoints);
1327 pointArray1.GetPoint(pIn1,npoints);
1328 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
1329 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
1330 DistortPoint(distPoint, sector);
1331 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
1332 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
1334 track.Rotate(alpha);
1335 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
1336 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
1337 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
1338 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
1339 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
1340 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
1341 pointArray0.AddPoint(npoints, &pIn0);
1342 pointArray1.AddPoint(npoints, &pIn1);
1344 if (npoints>=npoints0) break;
1346 if (npoints<npoints0/2) return 0;
1350 AliExternalTrackParam *track0=0;
1351 AliExternalTrackParam *track1=0;
1352 AliTrackPoint point1,point2,point3;
1353 if (dir==1) { //make seed inner
1354 pointArray0.GetPoint(point1,1);
1355 pointArray0.GetPoint(point2,30);
1356 pointArray0.GetPoint(point3,60);
1358 if (dir==-1){ //make seed outer
1359 pointArray0.GetPoint(point1,npoints-60);
1360 pointArray0.GetPoint(point2,npoints-30);
1361 pointArray0.GetPoint(point3,npoints-1);
1363 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
1364 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
1366 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
1367 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
1371 pointArray0.GetPoint(pIn0,ipoint);
1372 pointArray1.GetPoint(pIn1,ipoint);
1373 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
1374 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
1376 if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1377 if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1378 if (TMath::Abs(track0->GetZ())>kMaxZ) break;
1379 if (TMath::Abs(track0->GetX())>kMaxR) break;
1380 if (TMath::Abs(track1->GetZ())>kMaxZ) break;
1381 if (TMath::Abs(track1->GetX())>kMaxR) break;
1383 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
1385 Double_t pointPos[2]={0,0};
1386 Double_t pointCov[3]={0,0,0};
1387 pointPos[0]=prot0.GetY();//local y
1388 pointPos[1]=prot0.GetZ();//local z
1389 pointCov[0]=prot0.GetCov()[3];//simay^2
1390 pointCov[1]=prot0.GetCov()[4];//sigmayz
1391 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
1392 if (!track0->Update(pointPos,pointCov)) break;
1394 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
1395 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
1396 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
1398 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
1399 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
1400 pointCov[0]=prot1.GetCov()[3];//simay^2
1401 pointCov[1]=prot1.GetCov()[4];//sigmayz
1402 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
1403 if (!track1->Update(pointPos,pointCov)) break;
1407 if (npoints2<npoints) return 0;
1408 AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
1409 track1->Rotate(track0->GetAlpha());
1410 AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
1412 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
1413 "point0.="<<&pointArray0<< // points
1414 "point1.="<<&pointArray1<< // distorted points
1415 "trackIn.="<<&track<< // original track
1416 "track0.="<<track0<< // fitted track
1417 "track1.="<<track1<< // fitted distorted track
1419 new(&trackIn) AliExternalTrackParam(*track0);
1428 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
1430 // create the distortion tree on a mesh with granularity given by step
1431 // return the tree with distortions at given position
1432 // Map is created on the mesh with given step size
1434 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
1436 for (Double_t x= -250; x<250; x+=step){
1437 for (Double_t y= -250; y<250; y+=step){
1438 Double_t r = TMath::Sqrt(x*x+y*y);
1440 if (r>250) continue;
1441 for (Double_t z= -250; z<250; z+=step){
1442 Int_t roc=(z>0)?0:18;
1446 Double_t phi = TMath::ATan2(y,x);
1447 DistortPoint(xyz,roc);
1448 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1449 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
1450 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
1451 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
1452 Double_t dx = xyz[0]-x;
1453 Double_t dy = xyz[1]-y;
1454 Double_t dz = xyz[2]-z;
1456 Double_t drphi=(phi1-phi)*r;
1457 (*pcstream)<<"distortion"<<
1458 "x="<<x<< // original position
1463 "x1="<<xyz[0]<< // distorted position
1469 "dx="<<dx<< // delta position
1479 TFile f(Form("correction%s.root",GetName()));
1480 TTree * tree = (TTree*)f.Get("distortion");
1481 TTree * tree2= tree->CopyTree("1");
1482 tree2->SetName(Form("dist%s",GetName()));
1483 tree2->SetDirectory(0);
1491 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
1494 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1495 // calculates partial distortions
1496 // Partial distortion is stored in the resulting tree
1497 // Output is storred in the file distortion_<dettype>_<partype>.root
1498 // Partial distortion is stored with the name given by correction name
1501 // Parameters of function:
1502 // input - input tree
1503 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex
1504 // ppype - parameter type
1505 // corrArray - array with partial corrections
1506 // step - skipe entries - if 1 all entries processed - it is slow
1507 // debug 0 if debug on also space points dumped - it is slow
1509 const Double_t kMaxSnp = 0.85;
1510 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1511 // const Double_t kB2C=-0.299792458e-3;
1512 const Int_t kMinEntries=50;
1513 Double_t phi,theta, snp, mean,rms, entries;
1514 tinput->SetBranchAddress("theta",&theta);
1515 tinput->SetBranchAddress("phi", &phi);
1516 tinput->SetBranchAddress("snp",&snp);
1517 tinput->SetBranchAddress("mean",&mean);
1518 tinput->SetBranchAddress("rms",&rms);
1519 tinput->SetBranchAddress("entries",&entries);
1520 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
1522 Int_t nentries=tinput->GetEntries();
1523 Int_t ncorr=corrArray->GetEntries();
1524 Double_t corrections[100]={0}; //
1526 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1529 if (dtype==0) {refX=85.; dir=-1;}
1530 if (dtype==1) {refX=275.; dir=1;}
1531 if (dtype==2) {refX=85.; dir=-1;}
1532 if (dtype==3) {refX=360.; dir=-1;}
1534 for (Int_t ientry=0; ientry<nentries; ientry+=step){
1535 tinput->GetEntry(ientry);
1536 if (TMath::Abs(snp)>kMaxSnp) continue;
1541 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
1542 Double_t bz=AliTrackerBase::GetBz();
1543 if (refX>10. && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
1544 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
1545 AliExternalTrackParam track(refX,phi,tPar,cov);
1549 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1550 if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature
1551 (*pcstream)<<"fit"<<
1552 "bz="<<bz<< // magnetic filed used
1553 "dtype="<<dtype<< // detector match type
1554 "ptype="<<ptype<< // parameter type
1555 "theta="<<theta<< // theta
1556 "phi="<<phi<< // phi
1557 "snp="<<snp<< // snp
1558 "mean="<<mean<< // mean dist value
1559 "rms="<<rms<< // rms
1560 "gx="<<xyz[0]<< // global position at reference
1561 "gy="<<xyz[1]<< // global position at reference
1562 "gz="<<xyz[2]<< // global position at reference
1563 "dRrec="<<dRrec<< // delta Radius in reconstruction
1564 "id="<<id<< // track id
1565 "entries="<<entries;// number of entries in bin
1567 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1568 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1569 corrections[icorr]=0;
1570 if (entries>kMinEntries){
1571 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1572 AliExternalTrackParam *trackOut = 0;
1573 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1574 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1575 if (dtype==0) {refX=85.; dir=-1;}
1576 if (dtype==1) {refX=275.; dir=1;}
1577 if (dtype==2) {refX=0; dir=-1;}
1578 if (dtype==3) {refX=360.; dir=-1;}
1581 AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1582 trackOut->Rotate(trackIn.GetAlpha());
1583 trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1585 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1588 corrections[icorr]=0;
1590 if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature
1593 (*pcstream)<<"fit"<<
1594 Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
1595 Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
1596 // for points it is neccessary
1598 (*pcstream)<<"fit"<<"\n";
1605 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1607 // Make a laser fit tree for global minimization
1609 const Double_t cutErrY=0.1;
1610 const Double_t cutErrZ=0.1;
1611 const Double_t kEpsilon=0.00000001;
1616 AliTPCLaserTrack *ltr=0;
1617 AliTPCLaserTrack::LoadTracks();
1618 tree->SetBranchAddress("dY.",&vecdY);
1619 tree->SetBranchAddress("dZ.",&vecdZ);
1620 tree->SetBranchAddress("eY.",&veceY);
1621 tree->SetBranchAddress("eZ.",&veceZ);
1622 tree->SetBranchAddress("LTr.",<r);
1623 Int_t entries= tree->GetEntries();
1624 TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1625 Double_t bz=AliTrackerBase::GetBz();
1628 for (Int_t ientry=0; ientry<entries; ientry++){
1629 tree->GetEntry(ientry);
1630 if (!ltr->GetVecGX()){
1631 ltr->UpdatePoints();
1633 TVectorD * delta= (itype==0)? vecdY:vecdZ;
1634 TVectorD * err= (itype==0)? veceY:veceZ;
1636 for (Int_t irow=0; irow<159; irow++){
1637 Int_t nentries = 1000;
1638 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1639 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1641 Double_t phi =(*ltr->GetVecPhi())[irow];
1642 Double_t theta =ltr->GetTgl();
1643 Double_t mean=delta->GetMatrixArray()[irow];
1644 Double_t gx=0,gy=0,gz=0;
1645 Double_t snp = (*ltr->GetVecP2())[irow];
1646 Double_t rms = 0.1+err->GetMatrixArray()[irow];
1647 gx = (*ltr->GetVecGX())[irow];
1648 gy = (*ltr->GetVecGY())[irow];
1649 gz = (*ltr->GetVecGZ())[irow];
1650 Int_t bundle= ltr->GetBundle();
1653 // get delta R used in reconstruction
1654 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
1655 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1656 const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1657 Double_t xyz0[3]={gx,gy,gz};
1658 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1660 // old ExB correction
1662 if(recoParam&&recoParam->GetUseExBCorrection()) {
1663 Double_t xyz1[3]={gx,gy,gz};
1664 calib->GetExB()->Correct(xyz0,xyz1);
1665 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1668 if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1669 Float_t xyz1[3]={gx,gy,gz};
1670 Int_t sector=(gz>0)?0:18;
1671 correction->CorrectPoint(xyz1, sector);
1672 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1677 (*pcstream)<<"fit"<<
1678 "bz="<<bz<< // magnetic filed used
1679 "dtype="<<dtype<< // detector match type
1680 "ptype="<<itype<< // parameter type
1681 "theta="<<theta<< // theta
1682 "phi="<<phi<< // phi
1683 "snp="<<snp<< // snp
1684 "mean="<<mean<< // mean dist value
1685 "rms="<<rms<< // rms
1686 "gx="<<gx<< // global position
1687 "gy="<<gy<< // global position
1688 "gz="<<gz<< // global position
1689 "dRrec="<<dRrec<< // delta Radius in reconstruction
1690 "id="<<bundle<< //bundle
1691 "entries="<<nentries;// number of entries in bin
1694 Double_t ky = TMath::Tan(TMath::ASin(snp));
1695 Int_t ncorr = corrArray->GetEntries();
1696 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
1697 Double_t phi0 = TMath::ATan2(gy,gx);
1698 Double_t distortions[1000]={0};
1699 Double_t distortionsR[1000]={0};
1700 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1701 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1702 Float_t distPoint[3]={gx,gy,gz};
1703 Int_t sector= (gz>0)? 0:18;
1705 corr->DistortPoint(distPoint, sector);
1707 // Double_t value=distPoint[2]-gz;
1709 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1710 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1711 Double_t drphi= r0*(phi1-phi0);
1712 Double_t dr = r1-r0;
1713 distortions[icorr] = drphi-ky*dr;
1714 distortionsR[icorr] = dr;
1716 (*pcstream)<<"fit"<<
1717 Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
1718 Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
1720 (*pcstream)<<"fit"<<"\n";
1728 void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
1730 // make a distortion map out ou fthe residual histogram
1731 // Results are written to the debug streamer - pcstream
1733 // his0 - input (4D) residual histogram
1734 // pcstream - file to write the tree
1736 // marian.ivanov@cern.ch
1737 const Int_t kMinEntries=50;
1738 Int_t nbins1=his0->GetAxis(1)->GetNbins();
1739 Int_t first1=his0->GetAxis(1)->GetFirst();
1740 Int_t last1 =his0->GetAxis(1)->GetLast();
1742 Double_t bz=AliTrackerBase::GetBz();
1743 Int_t idim[4]={0,1,2,3};
1744 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
1746 his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1747 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1748 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
1749 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
1750 Int_t first3 = his1->GetAxis(3)->GetFirst();
1751 Int_t last3 = his1->GetAxis(3)->GetLast();
1754 for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){ // axis 3 - local angle
1755 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1756 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1758 his1->GetAxis(3)->SetRangeUser(-1,1);
1761 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
1762 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1763 Int_t first2 = his3->GetAxis(2)->GetFirst();
1764 Int_t last2 = his3->GetAxis(2)->GetLast();
1766 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1767 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1768 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1769 TH1 * hisDelta = his3->Projection(0);
1771 Double_t entries = hisDelta->GetEntries();
1772 Double_t mean=0, rms=0;
1773 if (entries>kMinEntries){
1774 mean = hisDelta->GetMean();
1775 rms = hisDelta->GetRMS();
1777 (*pcstream)<<hname<<
1783 "entries="<<entries<<
1788 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1800 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1802 // Store object in the OCDB
1803 // By default the object is stored in the current directory
1804 // default comment consit of user name and the date
1806 TString ocdbStorage="";
1807 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1808 AliCDBMetaData *metaData= new AliCDBMetaData();
1809 metaData->SetObjectClassName("AliTPCCorrection");
1810 metaData->SetResponsible("Marian Ivanov");
1811 metaData->SetBeamPeriod(1);
1812 metaData->SetAliRootVersion("05-25-01"); //root version
1813 TString userName=gSystem->GetFromPipe("echo $USER");
1814 TString date=gSystem->GetFromPipe("date");
1816 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1817 if (comment) metaData->SetComment(comment);
1819 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1820 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1821 gStorage->Put(this, (*id1), metaData);
1825 void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
1827 // Fast method to simulate the influence of the given distortion on the vertex reconstruction
1830 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
1831 if (!magF) AliError("Magneticd field - not initialized");
1832 Double_t bz = magF->SolenoidField(); //field in kGauss
1833 printf("bz: %lf\n",bz);
1834 AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
1836 TObjArray aTrk; // Original Track array of Aside
1837 TObjArray daTrk; // Distorted Track array of A side
1838 UShort_t *aId = new UShort_t[nTracks]; // A side Track ID
1841 UShort_t *cId = new UShort_t [nTracks];
1843 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1844 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
1845 fpt.SetParameters(7.24,0.120);
1847 for(Int_t nt=0; nt<nTracks; nt++){
1848 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
1849 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
1850 Double_t pt = fpt.GetRandom(); // momentum for f1
1851 // printf("phi %lf eta %lf pt %lf\n",phi,eta,pt);
1853 if(gRandom->Rndm() < 0.5){
1859 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1861 pxyz[0]=pt*TMath::Cos(phi);
1862 pxyz[1]=pt*TMath::Sin(phi);
1863 pxyz[2]=pt*TMath::Tan(theta);
1864 Double_t cv[21]={0};
1865 AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
1869 AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir, NULL);
1871 if (pcstream) (*pcstream)<<"track"<<
1877 if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
1881 Int_t nn=aTrk.GetEntriesFast();
1884 }else if(( eta<-0.07 )&&( eta>-etaCuts )){
1888 Int_t nn=cTrk.GetEntriesFast();
1893 }// end of track loop
1895 vertexer->SetTPCMode();
1896 vertexer->SetConstraintOff();
1898 aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));
1899 avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
1900 cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));
1901 cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
1902 if (pcstream) (*pcstream)<<"vertex"<<
1903 "x="<<orgVertex[0]<<
1904 "y="<<orgVertex[1]<<
1905 "z="<<orgVertex[2]<<
1906 "av.="<<&aV<< // distorted vertex A side
1907 "cv.="<<&cV<< // distroted vertex C side
1908 "avO.="<<&avOrg<< // original vertex A side
1915 void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
1917 // make correction available for visualization using
1918 // TFormula, TFX and TTree::Draw
1919 // important in order to check corrections and also compute dervied variables
1920 // e.g correction partial derivatives
1922 // NOTE - class is not owner of correction
1924 if (!fgVisualCorrection) fgVisualCorrection=new TObjArray;
1925 if (position>=fgVisualCorrection->GetEntriesFast()) fgVisualCorrection->Expand(position*2);
1926 fgVisualCorrection->AddAt(corr, position);
1931 Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t r, Double_t kZ, Int_t axisType, Int_t corrType){
1933 // calculate the correction at given position - check the geffCorr
1935 if (!fgVisualCorrection) return 0;
1936 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1937 if (!corr) return 0;
1938 Double_t phi=sector*TMath::Pi()/9.;
1939 Double_t gx = r*TMath::Cos(phi);
1940 Double_t gy = r*TMath::Sin(phi);
1942 Int_t nsector=(gz>0) ? 0:18;
1946 Float_t distPoint[3]={gx,gy,gz};
1947 corr->DistortPoint(distPoint, nsector);
1948 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1949 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1950 Double_t phi0=TMath::ATan2(gy,gx);
1951 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1952 if (axisType==0) return r1-r0;
1953 if (axisType==1) return (phi1-phi0)*r0;
1954 if (axisType==2) return distPoint[2]-gz;
1958 Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
1960 // return correction at given x,y,z
1962 if (!fgVisualCorrection) return 0;
1963 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1964 if (!corr) return 0;
1965 Double_t phi0= TMath::ATan2(gy,gx);
1966 Int_t nsector=(gz>0) ? 0:18;
1967 Float_t distPoint[3]={gx,gy,gz};
1968 corr->DistortPoint(distPoint, nsector);
1969 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1970 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1971 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1972 if (axisType==0) return r1-r0;
1973 if (axisType==1) return (phi1-phi0)*r0;
1974 if (axisType==2) return distPoint[2]-gz;