]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCCorrection.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update to speed up the calculation (and the accuracy) to calculate space charge effec...
[u/mrichter/AliRoot.git] / TPC / AliTPCCorrection.cxx
CommitLineData
0116859c 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16////////////////////////////////////////////////////////////////////////////////
17// //
18// AliTPCCorrection class //
19// //
20// This class provides a general framework to deal with space point //
21// distortions. An correction class which inherits from here is for example //
22// AliTPCExBBShape or AliTPCExBTwist //
23// //
24// General functions are (for example): //
25// CorrectPoint(x,roc) where x is the vector of inital positions in //
26// cartesian coordinates and roc represents the Read Out chamber number //
27// according to the offline naming convention. The vector x is overwritten //
28// with the corrected coordinates. //
29// //
30// An alternative usage would be CorrectPoint(x,roc,dx), which leaves the //
31// vector x untouched, put returns the distortions via the vector dx //
32// //
33// The class allows "effective Omega Tau" corrections to be shifted to the //
34// single distortion classes. //
35// //
36// Note: This class is normally used via the class AliTPCComposedCorrection //
37// //
38// date: 27/04/2010 //
39// Authors: Magnus Mager, Stefan Rossegger, Jim Thomas //
40////////////////////////////////////////////////////////////////////////////////
be67055b 41#include "Riostream.h"
0116859c 42
43#include <TH2F.h>
44#include <TMath.h>
45#include <TROOT.h>
cf5b0aa0 46#include <TTreeStream.h>
ffab0c37 47#include <TTree.h>
48#include <TFile.h>
e527a1b9 49#include <TTimeStamp.h>
ffab0c37 50#include <AliCDBStorage.h>
51#include <AliCDBId.h>
52#include <AliCDBMetaData.h>
c9cbd2f2 53#include "TVectorD.h"
54#include "AliTPCParamSR.h"
7f4cb119 55
c9cbd2f2 56#include "AliTPCCorrection.h"
57#include "AliLog.h"
1b923461 58
1b923461 59#include "AliExternalTrackParam.h"
60#include "AliTrackPointArray.h"
61#include "TDatabasePDG.h"
62#include "AliTrackerBase.h"
63#include "AliTPCROC.h"
64#include "THnSparse.h"
65
c9cbd2f2 66#include "AliTPCLaserTrack.h"
67#include "AliESDVertex.h"
68#include "AliVertexerTracks.h"
69#include "TDatabasePDG.h"
70#include "TF1.h"
7f4cb119 71#include "TRandom.h"
c9cbd2f2 72
73#include "TDatabasePDG.h"
74
7f4cb119 75#include "AliTPCTransform.h"
76#include "AliTPCcalibDB.h"
77#include "AliTPCExB.h"
cf5b0aa0 78
c9cbd2f2 79#include "AliTPCRecoParam.h"
1b923461 80
0116859c 81
cf5b0aa0 82ClassImp(AliTPCCorrection)
83
f1817479 84
85TObjArray *AliTPCCorrection::fgVisualCorrection=0;
86// instance of correction for visualization
87
88
0116859c 89// FIXME: the following values should come from the database
c9cbd2f2 90const Double_t AliTPCCorrection::fgkTPCZ0 = 249.7; // nominal gating grid position
2b68ab9c 91const 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)
c9cbd2f2 93const Double_t AliTPCCorrection::fgkOFCRadius= 254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
94const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
95const Double_t AliTPCCorrection::fgkCathodeV = -100000.0; // Cathode Voltage (volts)
96const Double_t AliTPCCorrection::fgkGG = -70.0; // Gating Grid voltage (volts)
0116859c 97
c9cbd2f2 98const Double_t AliTPCCorrection::fgkdvdE = 0.0024; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
0116859c 99
c9cbd2f2 100const Double_t AliTPCCorrection::fgkEM = -1.602176487e-19/9.10938215e-31; // charge/mass in [C/kg]
101const Double_t AliTPCCorrection::fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)]
c9cbd2f2 102
0116859c 103
104AliTPCCorrection::AliTPCCorrection()
c9cbd2f2 105 : TNamed("correction_unity","unity"),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 106{
107 //
108 // default constructor
109 //
f1817479 110 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
c9cbd2f2 111
35ae345f 112 InitLookUpfulcrums();
c9cbd2f2 113
0116859c 114}
115
116AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
c9cbd2f2 117: TNamed(name,title),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 118{
119 //
120 // default constructor, that set the name and title
121 //
f1817479 122 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
c9cbd2f2 123
35ae345f 124 InitLookUpfulcrums();
c9cbd2f2 125
0116859c 126}
127
128AliTPCCorrection::~AliTPCCorrection() {
129 //
130 // virtual destructor
131 //
132}
133
134void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
135 //
136 // Corrects the initial coordinates x (cartesian coordinates)
137 // according to the given effect (inherited classes)
138 // roc represents the TPC read out chamber (offline numbering convention)
139 //
140 Float_t dx[3];
141 GetCorrection(x,roc,dx);
142 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
143}
144
145void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
146 //
147 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
148 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
149 // roc represents the TPC read out chamber (offline numbering convention)
150 //
151 Float_t dx[3];
152 GetCorrection(x,roc,dx);
153 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
154}
155
156void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
157 //
158 // Distorts the initial coordinates x (cartesian coordinates)
159 // according to the given effect (inherited classes)
160 // roc represents the TPC read out chamber (offline numbering convention)
161 //
162 Float_t dx[3];
163 GetDistortion(x,roc,dx);
164 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
165}
166
167void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
168 //
169 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
170 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
171 // roc represents the TPC read out chamber (offline numbering convention)
172 //
173 Float_t dx[3];
174 GetDistortion(x,roc,dx);
175 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
176}
177
178void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
179 //
180 // This function delivers the correction values dx in respect to the inital coordinates x
181 // roc represents the TPC read out chamber (offline numbering convention)
182 // Note: The dx is overwritten by the inherited effectice class ...
183 //
184 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
185}
186
187void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
188 //
189 // This function delivers the distortion values dx in respect to the inital coordinates x
190 // roc represents the TPC read out chamber (offline numbering convention)
191 //
192 GetCorrection(x,roc,dx);
193 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
194}
195
196void AliTPCCorrection::Init() {
197 //
198 // Initialization funtion (not used at the moment)
199 //
200}
201
e527a1b9 202void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
203 //
204 // Update function
205 //
206}
207
0116859c 208void AliTPCCorrection::Print(Option_t* /*option*/) const {
209 //
210 // Print function to check which correction classes are used
211 // option=="d" prints details regarding the setted magnitude
212 // option=="a" prints the C0 and C1 coefficents for calibration purposes
213 //
214 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
215}
216
534fd34a 217void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
0116859c 218 //
219 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
220 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
221 // calibration run
222 //
534fd34a 223 fT1=t1;
224 fT2=t2;
225 //SetOmegaTauT1T2(omegaTau, t1, t2);
0116859c 226}
227
228TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
229 //
230 // Simple plot functionality.
231 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
232 // in respect to position z within the XY plane.
233 // The histogramm has nx times ny entries.
234 //
c9cbd2f2 235 AliTPCParam* tpcparam = new AliTPCParamSR;
236
0116859c 237 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
238 nx,-250.,250.,ny,-250.,250.);
239 Float_t x[3],dx[3];
240 x[2]=z;
241 Int_t roc=z>0.?0:18; // FIXME
242 for (Int_t iy=1;iy<=ny;++iy) {
243 x[1]=h->GetYaxis()->GetBinCenter(iy);
244 for (Int_t ix=1;ix<=nx;++ix) {
245 x[0]=h->GetXaxis()->GetBinCenter(ix);
246 GetCorrection(x,roc,dx);
247 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
c9cbd2f2 248 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
0116859c 249 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
250 h->SetBinContent(ix,iy,r1-r0);
251 }
252 else
253 h->SetBinContent(ix,iy,0.);
254 }
255 }
c9cbd2f2 256 delete tpcparam;
0116859c 257 return h;
258}
259
260TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
261 //
262 // Simple plot functionality.
263 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
264 // in respect to position z within the XY plane.
265 // The histogramm has nx times ny entries.
266 //
267
c9cbd2f2 268 AliTPCParam* tpcparam = new AliTPCParamSR;
269
0116859c 270 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
271 nx,-250.,250.,ny,-250.,250.);
272 Float_t x[3],dx[3];
273 x[2]=z;
274 Int_t roc=z>0.?0:18; // FIXME
275 for (Int_t iy=1;iy<=ny;++iy) {
276 x[1]=h->GetYaxis()->GetBinCenter(iy);
277 for (Int_t ix=1;ix<=nx;++ix) {
278 x[0]=h->GetXaxis()->GetBinCenter(ix);
279 GetCorrection(x,roc,dx);
280 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
c9cbd2f2 281 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
0116859c 282 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
283 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
284
285 Float_t dphi=phi1-phi0;
286 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
287 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
288
289 h->SetBinContent(ix,iy,r0*dphi);
290 }
291 else
292 h->SetBinContent(ix,iy,0.);
293 }
294 }
c9cbd2f2 295 delete tpcparam;
296 return h;
297}
298
299TH2F* AliTPCCorrection::CreateHistoDZinXY(Float_t z,Int_t nx,Int_t ny) {
300 //
301 // Simple plot functionality.
302 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
303 // in respect to position z within the XY plane.
304 // The histogramm has nx times ny entries.
305 //
306
307 AliTPCParam* tpcparam = new AliTPCParamSR;
308
309 TH2F *h=CreateTH2F("dz_xy",GetTitle(),"x [cm]","y [cm]","dz [cm]",
310 nx,-250.,250.,ny,-250.,250.);
311 Float_t x[3],dx[3];
312 x[2]=z;
313 Int_t roc=z>0.?0:18; // FIXME
314 for (Int_t iy=1;iy<=ny;++iy) {
315 x[1]=h->GetYaxis()->GetBinCenter(iy);
316 for (Int_t ix=1;ix<=nx;++ix) {
317 x[0]=h->GetXaxis()->GetBinCenter(ix);
318 GetCorrection(x,roc,dx);
319 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
320 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
321 h->SetBinContent(ix,iy,dx[2]);
322 }
323 else
324 h->SetBinContent(ix,iy,0.);
325 }
326 }
327 delete tpcparam;
0116859c 328 return h;
329}
330
331TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
332 //
333 // Simple plot functionality.
334 // Returns a 2d hisogram which represents the corrections in r direction (dr)
335 // in respect to angle phi within the ZR plane.
336 // The histogramm has nx times ny entries.
337 //
338 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
339 nz,-250.,250.,nr,85.,250.);
340 Float_t x[3],dx[3];
341 for (Int_t ir=1;ir<=nr;++ir) {
342 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
343 x[0]=radius*TMath::Cos(phi);
344 x[1]=radius*TMath::Sin(phi);
345 for (Int_t iz=1;iz<=nz;++iz) {
346 x[2]=h->GetXaxis()->GetBinCenter(iz);
347 Int_t roc=x[2]>0.?0:18; // FIXME
348 GetCorrection(x,roc,dx);
349 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
350 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
351 h->SetBinContent(iz,ir,r1-r0);
352 }
353 }
0116859c 354 return h;
355
356}
357
358TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
359 //
360 // Simple plot functionality.
361 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
362 // in respect to angle phi within the ZR plane.
363 // The histogramm has nx times ny entries.
364 //
365 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
366 nz,-250.,250.,nr,85.,250.);
367 Float_t x[3],dx[3];
368 for (Int_t iz=1;iz<=nz;++iz) {
369 x[2]=h->GetXaxis()->GetBinCenter(iz);
370 Int_t roc=x[2]>0.?0:18; // FIXME
371 for (Int_t ir=1;ir<=nr;++ir) {
372 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
373 x[0]=radius*TMath::Cos(phi);
374 x[1]=radius*TMath::Sin(phi);
375 GetCorrection(x,roc,dx);
376 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
377 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
378 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
379
380 Float_t dphi=phi1-phi0;
381 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
382 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
383
384 h->SetBinContent(iz,ir,r0*dphi);
385 }
386 }
387 return h;
388}
389
c9cbd2f2 390TH2F* AliTPCCorrection::CreateHistoDZinZR(Float_t phi,Int_t nz,Int_t nr) {
391 //
392 // Simple plot functionality.
393 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
394 // in respect to angle phi within the ZR plane.
395 // The histogramm has nx times ny entries.
396 //
397 TH2F *h=CreateTH2F("dz_zr",GetTitle(),"z [cm]","r [cm]","dz [cm]",
398 nz,-250.,250.,nr,85.,250.);
399 Float_t x[3],dx[3];
400 for (Int_t ir=1;ir<=nr;++ir) {
401 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
402 x[0]=radius*TMath::Cos(phi);
403 x[1]=radius*TMath::Sin(phi);
404 for (Int_t iz=1;iz<=nz;++iz) {
405 x[2]=h->GetXaxis()->GetBinCenter(iz);
406 Int_t roc=x[2]>0.?0:18; // FIXME
407 GetCorrection(x,roc,dx);
408 h->SetBinContent(iz,ir,dx[2]);
409 }
410 }
411 return h;
412
413}
414
415
0116859c 416TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
417 const char *xlabel,const char *ylabel,const char *zlabel,
418 Int_t nbinsx,Double_t xlow,Double_t xup,
419 Int_t nbinsy,Double_t ylow,Double_t yup) {
420 //
421 // Helper function to create a 2d histogramm of given size
422 //
423
424 TString hname=name;
425 Int_t i=0;
426 if (gDirectory) {
427 while (gDirectory->FindObject(hname.Data())) {
428 hname =name;
429 hname+="_";
430 hname+=i;
431 ++i;
432 }
433 }
434 TH2F *h=new TH2F(hname.Data(),title,
435 nbinsx,xlow,xup,
436 nbinsy,ylow,yup);
437 h->GetXaxis()->SetTitle(xlabel);
438 h->GetYaxis()->SetTitle(ylabel);
439 h->GetZaxis()->SetTitle(zlabel);
440 h->SetStats(0);
441 return h;
442}
443
0116859c 444// Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
445
446void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
b1f0a2a5 447 const Double_t er[kNZ][kNR], Double_t &erValue ) {
0116859c 448 //
449 // Interpolate table - 2D interpolation
450 //
25732bff 451 Double_t saveEr[5] = {0,0,0,0,0};
0116859c 452
453 Search( kNZ, fgkZList, z, fJLow ) ;
454 Search( kNR, fgkRList, r, fKLow ) ;
455 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
456 if ( fKLow < 0 ) fKLow = 0 ;
457 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
458 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
459
460 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
b1f0a2a5 461 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
0116859c 462 }
b1f0a2a5 463 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
0116859c 464
465}
466
c9cbd2f2 467void AliTPCCorrection::Interpolate3DEdistortion( const Int_t order, const Double_t r, const Float_t phi, const Double_t z,
468 const Double_t er[kNZ][kNPhi][kNR], const Double_t ephi[kNZ][kNPhi][kNR], const Double_t ez[kNZ][kNPhi][kNR],
469 Double_t &erValue, Double_t &ephiValue, Double_t &ezValue) {
470 //
471 // Interpolate table - 3D interpolation
472 //
473
25732bff 474 Double_t saveEr[5]= {0,0,0,0,0};
475 Double_t savedEr[5]= {0,0,0,0,0} ;
476
477 Double_t saveEphi[5]= {0,0,0,0,0};
478 Double_t savedEphi[5]= {0,0,0,0,0} ;
479
480 Double_t saveEz[5]= {0,0,0,0,0};
481 Double_t savedEz[5]= {0,0,0,0,0} ;
c9cbd2f2 482
483 Search( kNZ, fgkZList, z, fILow ) ;
484 Search( kNPhi, fgkPhiList, z, fJLow ) ;
485 Search( kNR, fgkRList, r, fKLow ) ;
486
487 if ( fILow < 0 ) fILow = 0 ; // check if out of range
488 if ( fJLow < 0 ) fJLow = 0 ;
489 if ( fKLow < 0 ) fKLow = 0 ;
490
491 if ( fILow + order >= kNZ - 1 ) fILow = kNZ - 1 - order ;
492 if ( fJLow + order >= kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
493 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
494
495 for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
496 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
497 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r ) ;
498 saveEphi[j-fJLow] = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
499 saveEz[j-fJLow] = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r ) ;
500 }
501 savedEr[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi ) ;
502 savedEphi[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
503 savedEz[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi ) ;
504 }
505 erValue = Interpolate( &fgkZList[fILow], savedEr, order, z ) ;
506 ephiValue = Interpolate( &fgkZList[fILow], savedEphi, order, z ) ;
507 ezValue = Interpolate( &fgkZList[fILow], savedEz, order, z ) ;
508
509}
510
511Double_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
512 const Int_t nx, const Int_t ny, const Double_t xv[], const Double_t yv[],
513 const TMatrixD &array ) {
514 //
515 // Interpolate table (TMatrix format) - 2D interpolation
516 //
517
518 static Int_t jlow = 0, klow = 0 ;
25732bff 519 Double_t saveArray[5] = {0,0,0,0,0} ;
c9cbd2f2 520
521 Search( nx, xv, x, jlow ) ;
522 Search( ny, yv, y, klow ) ;
523 if ( jlow < 0 ) jlow = 0 ; // check if out of range
524 if ( klow < 0 ) klow = 0 ;
525 if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
526 if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
527
528 for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
529 {
530 Double_t *ajkl = &((TMatrixD&)array)(j,klow);
531 saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
532 }
533
534 return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
535
536}
537
538Double_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x, const Double_t y, const Double_t z,
539 const Int_t nx, const Int_t ny, const Int_t nz,
540 const Double_t xv[], const Double_t yv[], const Double_t zv[],
541 TMatrixD **arrayofArrays ) {
542 //
543 // Interpolate table (TMatrix format) - 3D interpolation
544 //
545
546 static Int_t ilow = 0, jlow = 0, klow = 0 ;
25732bff 547 Double_t saveArray[5]= {0,0,0,0,0};
548 Double_t savedArray[5]= {0,0,0,0,0} ;
c9cbd2f2 549
550 Search( nx, xv, x, ilow ) ;
551 Search( ny, yv, y, jlow ) ;
552 Search( nz, zv, z, klow ) ;
553
554 if ( ilow < 0 ) ilow = 0 ; // check if out of range
555 if ( jlow < 0 ) jlow = 0 ;
556 if ( klow < 0 ) klow = 0 ;
557
558 if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
559 if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
560 if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
561
562 for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
563 {
564 TMatrixD &table = *arrayofArrays[k] ;
565 for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
566 {
567 saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
568 }
569 savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
570 }
571 return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
572
573}
574
575
0116859c 576Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
b1f0a2a5 577 const Int_t order, const Double_t x ) {
0116859c 578 //
579 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
580 //
581
582 Double_t y ;
583 if ( order == 2 ) { // Quadratic Interpolation = 2
584 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
585 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
586 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
587 } else { // Linear Interpolation = 1
588 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
589 }
590
591 return (y);
592
593}
594
595
b1f0a2a5 596void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
0116859c 597 //
598 // Search an ordered table by starting at the most recently used point
599 //
600
601 Long_t middle, high ;
602 Int_t ascend = 0, increment = 1 ;
603
604 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
605
606 if ( low < 0 || low > n-1 ) {
607 low = -1 ; high = n ;
608 } else { // Ordered Search phase
609 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
610 if ( low == n-1 ) return ;
611 high = low + 1 ;
612 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
613 low = high ;
614 increment *= 2 ;
615 high = low + increment ;
616 if ( high > n-1 ) { high = n ; break ; }
617 }
618 } else {
619 if ( low == 0 ) { low = -1 ; return ; }
620 high = low - 1 ;
621 while ( (Int_t)( x < xArray[low] ) == ascend ) {
622 high = low ;
623 increment *= 2 ;
624 if ( increment >= high ) { low = -1 ; break ; }
625 else low = high - increment ;
626 }
627 }
628 }
629
630 while ( (high-low) != 1 ) { // Binary Search Phase
631 middle = ( high + low ) / 2 ;
632 if ( (Int_t)( x >= xArray[middle] ) == ascend )
633 low = middle ;
634 else
635 high = middle ;
636 }
637
638 if ( x == xArray[n-1] ) low = n-2 ;
639 if ( x == xArray[0] ) low = 0 ;
640
641}
642
35ae345f 643void AliTPCCorrection::InitLookUpfulcrums() {
644 //
645 // Initialization of interpolation points - for main look up table
646 // (course grid in the middle, fine grid on the borders)
647 //
648
649 AliTPCROC * roc = AliTPCROC::Instance();
650 const Double_t rLow = TMath::Floor(roc->GetPadRowRadii(0,0))-1; // first padRow plus some margin
651
652 // fulcrums in R
653 fgkRList[0] = rLow;
654 for (Int_t i = 1; i<kNR; i++) {
655 fgkRList[i] = fgkRList[i-1] + 3.5; // 3.5 cm spacing
656 if (fgkRList[i]<90 ||fgkRList[i]>245)
657 fgkRList[i] = fgkRList[i-1] + 0.5; // 0.5 cm spacing
658 else if (fgkRList[i]<100 || fgkRList[i]>235)
659 fgkRList[i] = fgkRList[i-1] + 1.5; // 1.5 cm spacing
660 else if (fgkRList[i]<120 || fgkRList[i]>225)
661 fgkRList[i] = fgkRList[i-1] + 2.5; // 2.5 cm spacing
662 }
663
664 // fulcrums in Z
665 fgkZList[0] = -249.5;
666 fgkZList[kNZ-1] = 249.5;
667 for (Int_t j = 1; j<kNZ/2; j++) {
668 fgkZList[j] = fgkZList[j-1];
669 if (TMath::Abs(fgkZList[j])< 0.15)
670 fgkZList[j] = fgkZList[j-1] + 0.09; // 0.09 cm spacing
671 else if(TMath::Abs(fgkZList[j])< 0.6)
672 fgkZList[j] = fgkZList[j-1] + 0.4; // 0.4 cm spacing
673 else if (TMath::Abs(fgkZList[j])< 2.5 || TMath::Abs(fgkZList[j])>248)
674 fgkZList[j] = fgkZList[j-1] + 0.5; // 0.5 cm spacing
675 else if (TMath::Abs(fgkZList[j])<10 || TMath::Abs(fgkZList[j])>235)
676 fgkZList[j] = fgkZList[j-1] + 1.5; // 1.5 cm spacing
677 else if (TMath::Abs(fgkZList[j])<25 || TMath::Abs(fgkZList[j])>225)
678 fgkZList[j] = fgkZList[j-1] + 2.5; // 2.5 cm spacing
679 else
680 fgkZList[j] = fgkZList[j-1] + 4; // 4 cm spacing
681
682 fgkZList[kNZ-j-1] = -fgkZList[j];
683 }
684
685 // fulcrums in phi
686 for (Int_t k = 0; k<kNPhi; k++)
687 fgkPhiList[k] = TMath::TwoPi()*k/(kNPhi-1);
688
689
690}
691
692
c9cbd2f2 693void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
694 TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
695 const Int_t rows, const Int_t columns, const Int_t iterations,
696 const Bool_t rocDisplacement ) {
1b923461 697 //
698 // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
699 //
700 // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
701 // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
702 // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
703 // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
704 // you wish to solve Laplaces equation however it should not contain random numbers or you will get
705 // random numbers back as a solution.
706 // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
707 // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
708 // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
709 // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
710 // number of points and solves the problem again. This happens several times until the maximum number
711 // of points has been included in the array.
712 //
713 // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
714 // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
715 //
c9cbd2f2 716 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
717 //
1b923461 718 // Original code by Jim Thomas (STAR TPC Collaboration)
719 //
720
721 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
722
723 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
724 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
725 const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
726
727 TMatrixD arrayEr(rows,columns) ;
728 TMatrixD arrayEz(rows,columns) ;
729
730 //Check that number of rows and columns is suitable for a binary expansion
731
732 if ( !IsPowerOfTwo(rows-1) ) {
733 AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
734 return;
735 }
736 if ( !IsPowerOfTwo(columns-1) ) {
737 AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
738 return;
739 }
740
741 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
742 // Allow for different size grid spacing in R and Z directions
743 // Use a binary expansion of the size of the matrix to speed up the solution of the problem
744
745 Int_t iOne = (rows-1)/4 ;
746 Int_t jOne = (columns-1)/4 ;
747 // Solve for N in 2**N, add one.
748 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
749
750 for ( Int_t count = 0 ; count < loops ; count++ ) {
751 // Loop while the matrix expands & the resolution increases.
752
753 Float_t tempGridSizeR = gridSizeR * iOne ;
754 Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
755 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
756
757 // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
758 std::vector<float> coef1(rows) ;
759 std::vector<float> coef2(rows) ;
760
761 for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
762 Float_t radius = fgkIFCRadius + i*gridSizeR ;
763 coef1[i] = 1.0 + tempGridSizeR/(2*radius);
764 coef2[i] = 1.0 - tempGridSizeR/(2*radius);
765 }
766
767 TMatrixD sumChargeDensity(rows,columns) ;
768
769 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
770 Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
771 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
772 if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
773 else {
774 // Add up all enclosed charge density contributions within 1/2 unit in all directions
775 Float_t weight = 0.0 ;
776 Float_t sum = 0.0 ;
777 sumChargeDensity(i,j) = 0.0 ;
778 for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
779 for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
780 if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
781 else
782 weight = 1.0 ;
783 // Note that this is cylindrical geometry
784 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
785 sum += weight*radius ;
786 }
787 }
788 sumChargeDensity(i,j) /= sum ;
789 }
790 sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
791 }
792 }
793
794 for ( Int_t k = 1 ; k <= iterations; k++ ) {
795 // Solve Poisson's Equation
796 // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
797 // as interations increase.
798 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
799 Float_t overRelaxM1 = overRelax - 1.0 ;
800 Float_t overRelaxtempFourth, overRelaxcoef5 ;
801 overRelaxtempFourth = overRelax * tempFourth ;
802 overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
803
804 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
805 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
806
807 arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
808 + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
809 - overRelaxcoef5 * arrayV(i,j)
810 + coef1[i] * arrayV(i+iOne,j)
811 + sumChargeDensity(i,j)
812 ) * overRelaxtempFourth;
813 }
814 }
815
816 if ( k == iterations ) {
817 // After full solution is achieved, copy low resolution solution into higher res array
818 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
819 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
820
821 if ( iOne > 1 ) {
822 arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
823 if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
824 }
825 if ( jOne > 1 ) {
826 arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
827 if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
828 }
829 if ( iOne > 1 && jOne > 1 ) {
830 arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
831 if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
832 if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
833 // Note that this leaves a point at the upper left and lower right corners uninitialized.
834 // -> Not a big deal.
835 }
836
837 }
838 }
839 }
840
841 }
842
843 iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
844 jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
845
c9cbd2f2 846 sumChargeDensity.Clear();
1b923461 847 }
848
849 // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
850 for ( Int_t j = 0 ; j < columns ; j++ ) {
851 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
852 arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
853 arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
854 }
855
856 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
857 for ( Int_t i = 0 ; i < rows ; i++) {
858 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
859 arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
860 arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
861 }
862
863 for ( Int_t i = 0 ; i < rows ; i++) {
864 // Note: go back and compare to old version of this code. See notes below.
865 // JT Test ... attempt to divide by real Ez not Ez to first order
866 for ( Int_t j = 0 ; j < columns ; j++ ) {
867 arrayEz(i,j) += ezField;
868 // This adds back the overall Z gradient of the field (main E field component)
869 }
870 // Warning: (-=) assumes you are using an error potetial without the overall Field included
871 }
872
873 // Integrate Er/Ez from Z to zero
874 for ( Int_t j = 0 ; j < columns ; j++ ) {
875 for ( Int_t i = 0 ; i < rows ; i++ ) {
c9cbd2f2 876
1b923461 877 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
878 arrayErOverEz(i,j) = 0.0 ;
c9cbd2f2 879 arrayDeltaEz(i,j) = 0.0 ;
880
1b923461 881 for ( Int_t k = j ; k < columns ; k++ ) {
882 arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
c9cbd2f2 883 arrayDeltaEz(i,j) += index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
1b923461 884 if ( index != 4 ) index = 4; else index = 2 ;
885 }
c9cbd2f2 886 if ( index == 4 ) {
887 arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
888 arrayDeltaEz(i,j) -= (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
889 }
890 if ( index == 2 ) {
891 arrayErOverEz(i,j) += (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
892 -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
893 arrayDeltaEz(i,j) += (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
894 -2.5*(arrayEz(i,columns-1)-ezField));
895 }
896 if ( j == columns-2 ) {
897 arrayErOverEz(i,j) = (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
898 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
899 arrayDeltaEz(i,j) = (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
900 +1.5*(arrayEz(i,columns-1)-ezField) ) ;
901 }
902 if ( j == columns-1 ) {
903 arrayErOverEz(i,j) = 0.0 ;
904 arrayDeltaEz(i,j) = 0.0 ;
905 }
1b923461 906 }
907 }
908
c9cbd2f2 909 // calculate z distortion from the integrated Delta Ez residuals
910 // and include the aquivalence (Volt to cm) of the ROC shift !!
911
912 for ( Int_t j = 0 ; j < columns ; j++ ) {
913 for ( Int_t i = 0 ; i < rows ; i++ ) {
914
915 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
916 arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
917
918 // ROC Potential in cm aquivalent
919 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
920 if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift; // add the ROC misaligment
921
922 }
923 }
924
925 arrayEr.Clear();
926 arrayEz.Clear();
927
1b923461 928}
929
c9cbd2f2 930void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
931 TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
932 const Int_t rows, const Int_t columns, const Int_t phislices,
933 const Float_t deltaphi, const Int_t iterations, const Int_t symmetry,
934 Bool_t rocDisplacement ) {
935 //
936 // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
937 //
938 // NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
939 // The number of rows and COLUMNS can be different.
940 //
941 // ROWS == 2**M + 1
942 // COLUMNS == 2**N + 1
943 // PHISLICES == Arbitrary but greater than 3
944 //
945 // DeltaPhi in Radians
946 //
947 // SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
948 // = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
949 //
950 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
951
952 const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
953
954 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
955 const Float_t gridSizePhi = deltaphi ;
956 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
957 const Float_t ratioPhi = gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
958 const Float_t ratioZ = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
959
960 TMatrixD arrayE(rows,columns) ;
961
962 // Check that the number of rows and columns is suitable for a binary expansion
963 if ( !IsPowerOfTwo((rows-1)) ) {
964 AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
965 return; }
966 if ( !IsPowerOfTwo((columns-1)) ) {
967 AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
968 return; }
969 if ( phislices <= 3 ) {
970 AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
971 return; }
972 if ( phislices > 1000 ) {
973 AliError("Poisson3D phislices > 1000 is not allowed (nor wise) ");
974 return; }
975
976 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
977 // Allow for different size grid spacing in R and Z directions
978 // Use a binary expansion of the matrix to speed up the solution of the problem
979
980 Int_t loops, mplus, mminus, signplus, signminus ;
981 Int_t ione = (rows-1)/4 ;
982 Int_t jone = (columns-1)/4 ;
983 loops = TMath::Max(ione, jone) ; // Calculate the number of loops for the binary expansion
984 loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ; // Solve for N in 2**N
985
986 TMatrixD* arrayofSumChargeDensities[1000] ; // Create temporary arrays to store low resolution charge arrays
987
988 for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
989
990 for ( Int_t count = 0 ; count < loops ; count++ ) { // START the master loop and do the binary expansion
991
992 Float_t tempgridSizeR = gridSizeR * ione ;
993 Float_t tempratioPhi = ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
994 Float_t tempratioZ = ratioZ * ione * ione / ( jone * jone ) ;
995
996 std::vector<float> coef1(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
997 std::vector<float> coef2(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
998 std::vector<float> coef3(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
999 std::vector<float> coef4(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1000
1001 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1002 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1003 coef1[i] = 1.0 + tempgridSizeR/(2*radius);
1004 coef2[i] = 1.0 - tempgridSizeR/(2*radius);
1005 coef3[i] = tempratioPhi/(radius*radius);
1006 coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
1007 }
1008
1009 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1010 TMatrixD &chargeDensity = *arrayofChargeDensities[m] ;
1011 TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
1012 for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
1013 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1014 for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
1015 if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1016 else { // Add up all enclosed charge density contributions within 1/2 unit in all directions
1017 Float_t weight = 0.0 ;
1018 Float_t sum = 0.0 ;
1019 sumChargeDensity(i,j) = 0.0 ;
1020 for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
1021 for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
1022 if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
1023 else
1024 weight = 1.0 ;
1025 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1026 sum += weight*radius ;
1027 }
1028 }
1029 sumChargeDensity(i,j) /= sum ;
1030 }
1031 sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
1032 }
1033 }
1034 }
1035
1036 for ( Int_t k = 1 ; k <= iterations; k++ ) {
1037
1038 // over-relaxation index, >= 1 but < 2
1039 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1040 Float_t overRelaxM1 = overRelax - 1.0 ;
1041
1042 std::vector<float> overRelaxcoef4(rows) ; // Do this the standard C++ way to avoid gcc extensions
1043 std::vector<float> overRelaxcoef5(rows) ; // Do this the standard C++ way to avoid gcc extensions
1044
1045 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1046 overRelaxcoef4[i] = overRelax * coef4[i] ;
1047 overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
1048 }
1049
1050 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1051
1052 mplus = m + 1; signplus = 1 ;
1053 mminus = m - 1 ; signminus = 1 ;
1054 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1055 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1056 if ( mminus < 0 ) mminus = 1 ;
1057 }
1058 else if (symmetry==-1) { // Anti-symmetry in phi
1059 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1060 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1061 }
1062 else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1063 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1064 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1065 }
1066 TMatrixD& arrayV = *arrayofArrayV[m] ;
1067 TMatrixD& arrayVP = *arrayofArrayV[mplus] ;
1068 TMatrixD& arrayVM = *arrayofArrayV[mminus] ;
1069 TMatrixD& sumChargeDensity = *arrayofSumChargeDensities[m] ;
1070
1071 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1072 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1073
1074 arrayV(i,j) = ( coef2[i] * arrayV(i-ione,j)
1075 + tempratioZ * ( arrayV(i,j-jone) + arrayV(i,j+jone) )
1076 - overRelaxcoef5[i] * arrayV(i,j)
1077 + coef1[i] * arrayV(i+ione,j)
1078 + coef3[i] * ( signplus*arrayVP(i,j) + signminus*arrayVM(i,j) )
1079 + sumChargeDensity(i,j)
1080 ) * overRelaxcoef4[i] ;
1081 // Note: over-relax the solution at each step. This speeds up the convergance.
1082
1083 }
1084 }
1085
1086 if ( k == iterations ) { // After full solution is achieved, copy low resolution solution into higher res array
1087 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1088 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1089
1090 if ( ione > 1 ) {
1091 arrayV(i+ione/2,j) = ( arrayV(i+ione,j) + arrayV(i,j) ) / 2 ;
1092 if ( i == ione ) arrayV(i-ione/2,j) = ( arrayV(0,j) + arrayV(ione,j) ) / 2 ;
1093 }
1094 if ( jone > 1 ) {
1095 arrayV(i,j+jone/2) = ( arrayV(i,j+jone) + arrayV(i,j) ) / 2 ;
1096 if ( j == jone ) arrayV(i,j-jone/2) = ( arrayV(i,0) + arrayV(i,jone) ) / 2 ;
1097 }
1098 if ( ione > 1 && jone > 1 ) {
1099 arrayV(i+ione/2,j+jone/2) = ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
1100 if ( i == ione ) arrayV(i-ione/2,j-jone/2) = ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
1101 if ( j == jone ) arrayV(i-ione/2,j-jone/2) = ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
1102 // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
1103 }
1104 }
1105 }
1106 }
1107
1108 }
1109 }
1110
1111 ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
1112 jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
1113
1114 }
1115
1116 //Differentiate V(r) and solve for E(r) using special equations for the first and last row
1117 //Integrate E(r)/E(z) from point of origin to pad plane
1118
1119 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1120 TMatrixD& arrayV = *arrayofArrayV[m] ;
1121 TMatrixD& eroverEz = *arrayofEroverEz[m] ;
1122
1123 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1124
1125 // Differentiate in R
1126 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1127 arrayE(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1128 arrayE(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1129 // Integrate over Z
1130 for ( Int_t i = 0 ; i < rows ; i++ ) {
1131 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1132 eroverEz(i,j) = 0.0 ;
1133 for ( Int_t k = j ; k < columns ; k++ ) {
1134
1135 eroverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1136 if ( index != 4 ) index = 4; else index = 2 ;
1137 }
1138 if ( index == 4 ) eroverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1139 if ( index == 2 ) eroverEz(i,j) +=
1140 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1141 if ( j == columns-2 ) eroverEz(i,j) =
1142 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1143 if ( j == columns-1 ) eroverEz(i,j) = 0.0 ;
1144 }
1145 }
1146 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1147 // eroverEz.Draw("surf") ; } // JT test
1148 }
1149
1150 //Differentiate V(r) and solve for E(phi)
1151 //Integrate E(phi)/E(z) from point of origin to pad plane
1152
1153 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1154
1155 mplus = m + 1; signplus = 1 ;
1156 mminus = m - 1 ; signminus = 1 ;
1157 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1158 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1159 if ( mminus < 0 ) mminus = 1 ;
1160 }
1161 else if (symmetry==-1) { // Anti-symmetry in phi
1162 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1163 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1164 }
1165 else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
1166 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1167 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1168 }
1169 TMatrixD &arrayVP = *arrayofArrayV[mplus] ;
1170 TMatrixD &arrayVM = *arrayofArrayV[mminus] ;
1171 TMatrixD &ePhioverEz = *arrayofEPhioverEz[m] ;
1172 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1173 // Differentiate in Phi
1174 for ( Int_t i = 0 ; i < rows ; i++ ) {
1175 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1176 arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
1177 }
1178 // Integrate over Z
1179 for ( Int_t i = 0 ; i < rows ; i++ ) {
1180 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1181 ePhioverEz(i,j) = 0.0 ;
1182 for ( Int_t k = j ; k < columns ; k++ ) {
1183
1184 ePhioverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1185 if ( index != 4 ) index = 4; else index = 2 ;
1186 }
1187 if ( index == 4 ) ePhioverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1188 if ( index == 2 ) ePhioverEz(i,j) +=
1189 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1190 if ( j == columns-2 ) ePhioverEz(i,j) =
1191 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1192 if ( j == columns-1 ) ePhioverEz(i,j) = 0.0 ;
1193 }
1194 }
1195 // if ( m == 5 ) { TCanvas* c2 = new TCanvas("arrayE","arrayE",50,50,840,600) ; c2 -> cd() ;
1196 // arrayE.Draw("surf") ; } // JT test
1197 }
1198
1199
1200 // Differentiate V(r) and solve for E(z) using special equations for the first and last row
1201 // Integrate (E(z)-Ezstd) from point of origin to pad plane
1202
1203 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1204 TMatrixD& arrayV = *arrayofArrayV[m] ;
1205 TMatrixD& deltaEz = *arrayofDeltaEz[m] ;
1206
1207 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1208 for ( Int_t i = 0 ; i < rows ; i++) {
1209 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1210 arrayE(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1211 arrayE(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1212 }
1213
1214 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1215 // Integrate over Z
1216 for ( Int_t i = 0 ; i < rows ; i++ ) {
1217 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1218 deltaEz(i,j) = 0.0 ;
1219 for ( Int_t k = j ; k < columns ; k++ ) {
1220 deltaEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k) ;
1221 if ( index != 4 ) index = 4; else index = 2 ;
1222 }
1223 if ( index == 4 ) deltaEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1) ;
1224 if ( index == 2 ) deltaEz(i,j) +=
1225 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
1226 if ( j == columns-2 ) deltaEz(i,j) =
1227 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
1228 if ( j == columns-1 ) deltaEz(i,j) = 0.0 ;
1229 }
1230 }
1231 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1232 // eroverEz.Draw("surf") ; } // JT test
1233
1234 // calculate z distortion from the integrated Delta Ez residuals
1235 // and include the aquivalence (Volt to cm) of the ROC shift !!
1236
1237 for ( Int_t j = 0 ; j < columns ; j++ ) {
1238 for ( Int_t i = 0 ; i < rows ; i++ ) {
1239
1240 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1241 deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
1242
1243 // ROC Potential in cm aquivalent
1244 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
1245 if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift; // add the ROC misaligment
1246
1247 }
1248 }
1249
1250 } // end loop over phi
1251
1252
1253
1254 for ( Int_t k = 0 ; k < phislices ; k++ )
1255 {
1256 arrayofSumChargeDensities[k]->Delete() ;
1257 }
1258
1259
1260
1261 arrayE.Clear();
1262}
1b923461 1263
1264
710bda39 1265Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1b923461 1266 //
1267 // Helperfunction: Check if integer is a power of 2
1268 //
1269 Int_t j = 0;
1270 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
1271 if ( j == 1 ) return(1) ; // True
1272 return(0) ; // False
1273}
1274
cf5b0aa0 1275
b1f0a2a5 1276AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
cf5b0aa0 1277 //
1278 // Fit the track parameters - without and with distortion
1279 // 1. Space points in the TPC are simulated along the trajectory
1280 // 2. Space points distorted
1281 // 3. Fits the non distorted and distroted track to the reference plane at refX
1282 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
1283 //
1284 // trackIn - input track parameters
1285 // refX - reference X to fit the track
1286 // dir - direction - out=1 or in=-1
1287 // pcstream - debug streamer to check the results
1288 //
cad404e1 1289 // see AliExternalTrackParam.h documentation:
1290 // track1.fP[0] - local y (rphi)
1291 // track1.fP[1] - z
1292 // track1.fP[2] - sinus of local inclination angle
1293 // track1.fP[3] - tangent of deep angle
1294 // track1.fP[4] - 1/pt
1b923461 1295
cf5b0aa0 1296 AliTPCROC * roc = AliTPCROC::Instance();
1297 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
1298 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
1299 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
cf5b0aa0 1300 const Double_t kMaxSnp = 0.85;
1301 const Double_t kSigmaY=0.1;
1302 const Double_t kSigmaZ=0.1;
ca58ed4e 1303 const Double_t kMaxR=500;
1304 const Double_t kMaxZ=500;
cf5b0aa0 1305 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
ca58ed4e 1306 Int_t npoints1=0;
1307 Int_t npoints2=0;
cf5b0aa0 1308
be67055b 1309 AliExternalTrackParam track(trackIn); //
cf5b0aa0 1310 // generate points
1311 AliTrackPointArray pointArray0(npoints0);
1312 AliTrackPointArray pointArray1(npoints0);
1313 Double_t xyz[3];
ca58ed4e 1314 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 1315 //
1316 // simulate the track
1317 Int_t npoints=0;
1318 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
1319 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
ca58ed4e 1320 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 1321 track.GetXYZ(xyz);
ca58ed4e 1322 xyz[0]+=gRandom->Gaus(0,0.00005);
1323 xyz[1]+=gRandom->Gaus(0,0.00005);
1324 xyz[2]+=gRandom->Gaus(0,0.00005);
1325 if (TMath::Abs(track.GetZ())>kMaxZ) break;
1326 if (TMath::Abs(track.GetX())>kMaxR) break;
cf5b0aa0 1327 AliTrackPoint pIn0; // space point
1328 AliTrackPoint pIn1;
ffab0c37 1329 Int_t sector= (xyz[2]>0)? 0:18;
cf5b0aa0 1330 pointArray0.GetPoint(pIn0,npoints);
1331 pointArray1.GetPoint(pIn1,npoints);
1332 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
1333 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
ffab0c37 1334 DistortPoint(distPoint, sector);
cf5b0aa0 1335 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
1336 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
1337 //
1338 track.Rotate(alpha);
1339 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
1340 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
1341 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
1342 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
1343 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
1344 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
1345 pointArray0.AddPoint(npoints, &pIn0);
1346 pointArray1.AddPoint(npoints, &pIn1);
1347 npoints++;
1348 if (npoints>=npoints0) break;
1349 }
7f4cb119 1350 if (npoints<npoints0/2) return 0;
cf5b0aa0 1351 //
1352 // refit track
1353 //
1354 AliExternalTrackParam *track0=0;
1355 AliExternalTrackParam *track1=0;
1356 AliTrackPoint point1,point2,point3;
1357 if (dir==1) { //make seed inner
1358 pointArray0.GetPoint(point1,1);
4486a91f 1359 pointArray0.GetPoint(point2,30);
1360 pointArray0.GetPoint(point3,60);
cf5b0aa0 1361 }
1362 if (dir==-1){ //make seed outer
4486a91f 1363 pointArray0.GetPoint(point1,npoints-60);
1364 pointArray0.GetPoint(point2,npoints-30);
cf5b0aa0 1365 pointArray0.GetPoint(point3,npoints-1);
1366 }
1367 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
1368 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
1369
cf5b0aa0 1370 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
8b63d99c 1371 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
cf5b0aa0 1372 //
1373 AliTrackPoint pIn0;
1374 AliTrackPoint pIn1;
1375 pointArray0.GetPoint(pIn0,ipoint);
1376 pointArray1.GetPoint(pIn1,ipoint);
1377 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
1378 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
1379 //
ca58ed4e 1380 if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1381 if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1382 if (TMath::Abs(track0->GetZ())>kMaxZ) break;
1383 if (TMath::Abs(track0->GetX())>kMaxR) break;
1384 if (TMath::Abs(track1->GetZ())>kMaxZ) break;
1385 if (TMath::Abs(track1->GetX())>kMaxR) break;
1386
8b63d99c 1387 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
cf5b0aa0 1388 //
1389 Double_t pointPos[2]={0,0};
1390 Double_t pointCov[3]={0,0,0};
1391 pointPos[0]=prot0.GetY();//local y
1392 pointPos[1]=prot0.GetZ();//local z
1393 pointCov[0]=prot0.GetCov()[3];//simay^2
1394 pointCov[1]=prot0.GetCov()[4];//sigmayz
1395 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
ca58ed4e 1396 if (!track0->Update(pointPos,pointCov)) break;
cf5b0aa0 1397 //
8b63d99c 1398 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
1399 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
1400 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
1401
0b736a46 1402 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
1403 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
cf5b0aa0 1404 pointCov[0]=prot1.GetCov()[3];//simay^2
1405 pointCov[1]=prot1.GetCov()[4];//sigmayz
1406 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
ca58ed4e 1407 if (!track1->Update(pointPos,pointCov)) break;
1408 npoints1++;
1409 npoints2++;
cf5b0aa0 1410 }
ca58ed4e 1411 if (npoints2<npoints) return 0;
8b63d99c 1412 AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
cf5b0aa0 1413 track1->Rotate(track0->GetAlpha());
ca58ed4e 1414 AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
cf5b0aa0 1415
cad404e1 1416 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
cf5b0aa0 1417 "point0.="<<&pointArray0<< // points
1418 "point1.="<<&pointArray1<< // distorted points
1419 "trackIn.="<<&track<< // original track
1420 "track0.="<<track0<< // fitted track
1421 "track1.="<<track1<< // fitted distorted track
1422 "\n";
be67055b 1423 new(&trackIn) AliExternalTrackParam(*track0);
cf5b0aa0 1424 delete track0;
1425 return track1;
1426}
1427
1428
ffab0c37 1429
1430
1431
1432TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
1433 //
1434 // create the distortion tree on a mesh with granularity given by step
1435 // return the tree with distortions at given position
1436 // Map is created on the mesh with given step size
1437 //
1438 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
1439 Float_t xyz[3];
1440 for (Double_t x= -250; x<250; x+=step){
1441 for (Double_t y= -250; y<250; y+=step){
1442 Double_t r = TMath::Sqrt(x*x+y*y);
1443 if (r<80) continue;
1444 if (r>250) continue;
1445 for (Double_t z= -250; z<250; z+=step){
1446 Int_t roc=(z>0)?0:18;
1447 xyz[0]=x;
1448 xyz[1]=y;
1449 xyz[2]=z;
1450 Double_t phi = TMath::ATan2(y,x);
1451 DistortPoint(xyz,roc);
1452 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1453 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
1454 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
1455 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
1456 Double_t dx = xyz[0]-x;
1457 Double_t dy = xyz[1]-y;
1458 Double_t dz = xyz[2]-z;
1459 Double_t dr=r1-r;
1460 Double_t drphi=(phi1-phi)*r;
1461 (*pcstream)<<"distortion"<<
1462 "x="<<x<< // original position
1463 "y="<<y<<
1464 "z="<<z<<
1465 "r="<<r<<
1466 "phi="<<phi<<
1467 "x1="<<xyz[0]<< // distorted position
1468 "y1="<<xyz[1]<<
1469 "z1="<<xyz[2]<<
1470 "r1="<<r1<<
1471 "phi1="<<phi1<<
1472 //
1473 "dx="<<dx<< // delta position
1474 "dy="<<dy<<
1475 "dz="<<dz<<
1476 "dr="<<dr<<
1477 "drphi="<<drphi<<
1478 "\n";
1479 }
1480 }
1481 }
1482 delete pcstream;
1483 TFile f(Form("correction%s.root",GetName()));
1484 TTree * tree = (TTree*)f.Get("distortion");
1485 TTree * tree2= tree->CopyTree("1");
1486 tree2->SetName(Form("dist%s",GetName()));
1487 tree2->SetDirectory(0);
1488 delete tree;
1489 return tree2;
1490}
1491
1492
1493
be67055b 1494
b1f0a2a5 1495void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
be67055b 1496 //
1497 // Make a fit tree:
1498 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1499 // calculates partial distortions
1500 // Partial distortion is stored in the resulting tree
1501 // Output is storred in the file distortion_<dettype>_<partype>.root
1502 // Partial distortion is stored with the name given by correction name
1503 //
1504 //
1505 // Parameters of function:
1506 // input - input tree
1507 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex
1508 // ppype - parameter type
1509 // corrArray - array with partial corrections
1510 // step - skipe entries - if 1 all entries processed - it is slow
1511 // debug 0 if debug on also space points dumped - it is slow
c9cbd2f2 1512
b322e06a 1513 const Double_t kMaxSnp = 0.85;
1514 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1515 // const Double_t kB2C=-0.299792458e-3;
7f4cb119 1516 const Int_t kMinEntries=50;
be67055b 1517 Double_t phi,theta, snp, mean,rms, entries;
1518 tinput->SetBranchAddress("theta",&theta);
1519 tinput->SetBranchAddress("phi", &phi);
1520 tinput->SetBranchAddress("snp",&snp);
1521 tinput->SetBranchAddress("mean",&mean);
1522 tinput->SetBranchAddress("rms",&rms);
1523 tinput->SetBranchAddress("entries",&entries);
1524 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
1525 //
1526 Int_t nentries=tinput->GetEntries();
1527 Int_t ncorr=corrArray->GetEntries();
7f4cb119 1528 Double_t corrections[100]={0}; //
be67055b 1529 Double_t tPar[5];
1530 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1531 Double_t refX=0;
1532 Int_t dir=0;
7f4cb119 1533 if (dtype==0) {refX=85.; dir=-1;}
1534 if (dtype==1) {refX=275.; dir=1;}
1535 if (dtype==2) {refX=85.; dir=-1;}
1536 if (dtype==3) {refX=360.; dir=-1;}
be67055b 1537 //
1538 for (Int_t ientry=0; ientry<nentries; ientry+=step){
1539 tinput->GetEntry(ientry);
7f4cb119 1540 if (TMath::Abs(snp)>kMaxSnp) continue;
be67055b 1541 tPar[0]=0;
1542 tPar[1]=theta*refX;
1543 tPar[2]=snp;
1544 tPar[3]=theta;
4486a91f 1545 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
8b63d99c 1546 Double_t bz=AliTrackerBase::GetBz();
4486a91f 1547 if (refX>10. && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
1548 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
7f4cb119 1549 AliExternalTrackParam track(refX,phi,tPar,cov);
1550 Double_t xyz[3];
1551 track.GetXYZ(xyz);
1552 Int_t id=0;
1553 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
0b736a46 1554 if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature
be67055b 1555 (*pcstream)<<"fit"<<
8b63d99c 1556 "bz="<<bz<< // magnetic filed used
be67055b 1557 "dtype="<<dtype<< // detector match type
1558 "ptype="<<ptype<< // parameter type
1559 "theta="<<theta<< // theta
1560 "phi="<<phi<< // phi
1561 "snp="<<snp<< // snp
1562 "mean="<<mean<< // mean dist value
1563 "rms="<<rms<< // rms
7f4cb119 1564 "gx="<<xyz[0]<< // global position at reference
1565 "gy="<<xyz[1]<< // global position at reference
1566 "gz="<<xyz[2]<< // global position at reference
1567 "dRrec="<<dRrec<< // delta Radius in reconstruction
1568 "id="<<id<< // track id
be67055b 1569 "entries="<<entries;// number of entries in bin
1570 //
1571 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1572 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1573 corrections[icorr]=0;
1574 if (entries>kMinEntries){
1575 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1576 AliExternalTrackParam *trackOut = 0;
1577 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1578 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
7f4cb119 1579 if (dtype==0) {refX=85.; dir=-1;}
1580 if (dtype==1) {refX=275.; dir=1;}
b1f0a2a5 1581 if (dtype==2) {refX=0; dir=-1;}
7f4cb119 1582 if (dtype==3) {refX=360.; dir=-1;}
b1f0a2a5 1583 //
7f4cb119 1584 if (trackOut){
1585 AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1586 trackOut->Rotate(trackIn.GetAlpha());
1587 trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1588 //
1589 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1590 delete trackOut;
1591 }else{
1592 corrections[icorr]=0;
1593 }
0b736a46 1594 if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature
be67055b 1595 }
7f4cb119 1596 Double_t dRdummy=0;
be67055b 1597 (*pcstream)<<"fit"<<
7f4cb119 1598 Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
1599 Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
1600 // for points it is neccessary
be67055b 1601 }
1602 (*pcstream)<<"fit"<<"\n";
1603 }
1604 delete pcstream;
1605}
1606
1607
1608
7f4cb119 1609void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1610 //
1611 // Make a laser fit tree for global minimization
1612 //
1613 const Double_t cutErrY=0.1;
1614 const Double_t cutErrZ=0.1;
1615 const Double_t kEpsilon=0.00000001;
1616 TVectorD *vecdY=0;
1617 TVectorD *vecdZ=0;
1618 TVectorD *veceY=0;
1619 TVectorD *veceZ=0;
1620 AliTPCLaserTrack *ltr=0;
1621 AliTPCLaserTrack::LoadTracks();
1622 tree->SetBranchAddress("dY.",&vecdY);
1623 tree->SetBranchAddress("dZ.",&vecdZ);
1624 tree->SetBranchAddress("eY.",&veceY);
1625 tree->SetBranchAddress("eZ.",&veceZ);
1626 tree->SetBranchAddress("LTr.",&ltr);
1627 Int_t entries= tree->GetEntries();
1628 TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1629 Double_t bz=AliTrackerBase::GetBz();
1630 //
1631
1632 for (Int_t ientry=0; ientry<entries; ientry++){
1633 tree->GetEntry(ientry);
1634 if (!ltr->GetVecGX()){
1635 ltr->UpdatePoints();
1636 }
1637 TVectorD * delta= (itype==0)? vecdY:vecdZ;
1638 TVectorD * err= (itype==0)? veceY:veceZ;
1639
1640 for (Int_t irow=0; irow<159; irow++){
1641 Int_t nentries = 1000;
1642 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1643 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1644 Int_t dtype=4;
1645 Double_t phi =(*ltr->GetVecPhi())[irow];
1646 Double_t theta =ltr->GetTgl();
1647 Double_t mean=delta->GetMatrixArray()[irow];
1648 Double_t gx=0,gy=0,gz=0;
1649 Double_t snp = (*ltr->GetVecP2())[irow];
1650 Double_t rms = 0.1+err->GetMatrixArray()[irow];
1651 gx = (*ltr->GetVecGX())[irow];
1652 gy = (*ltr->GetVecGY())[irow];
1653 gz = (*ltr->GetVecGZ())[irow];
1654 Int_t bundle= ltr->GetBundle();
1655 Double_t dRrec=0;
1656 //
1657 // get delta R used in reconstruction
1658 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
1659 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1660 const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1661 Double_t xyz0[3]={gx,gy,gz};
1662 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1663 //
1664 // old ExB correction
1665 //
1666 if(recoParam&&recoParam->GetUseExBCorrection()) {
1667 Double_t xyz1[3]={gx,gy,gz};
1668 calib->GetExB()->Correct(xyz0,xyz1);
1669 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1670 dRrec=oldR-newR;
1671 }
1672 if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1673 Float_t xyz1[3]={gx,gy,gz};
1674 Int_t sector=(gz>0)?0:18;
1675 correction->CorrectPoint(xyz1, sector);
1676 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1677 dRrec=oldR-newR;
1678 }
1679
1680
1681 (*pcstream)<<"fit"<<
1682 "bz="<<bz<< // magnetic filed used
1683 "dtype="<<dtype<< // detector match type
1684 "ptype="<<itype<< // parameter type
1685 "theta="<<theta<< // theta
1686 "phi="<<phi<< // phi
1687 "snp="<<snp<< // snp
1688 "mean="<<mean<< // mean dist value
1689 "rms="<<rms<< // rms
1690 "gx="<<gx<< // global position
1691 "gy="<<gy<< // global position
1692 "gz="<<gz<< // global position
1693 "dRrec="<<dRrec<< // delta Radius in reconstruction
1694 "id="<<bundle<< //bundle
1695 "entries="<<nentries;// number of entries in bin
1696 //
1697 //
1698 Double_t ky = TMath::Tan(TMath::ASin(snp));
1699 Int_t ncorr = corrArray->GetEntries();
1700 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
1701 Double_t phi0 = TMath::ATan2(gy,gx);
1702 Double_t distortions[1000]={0};
1703 Double_t distortionsR[1000]={0};
1704 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1705 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1706 Float_t distPoint[3]={gx,gy,gz};
1707 Int_t sector= (gz>0)? 0:18;
1708 if (r0>80){
1709 corr->DistortPoint(distPoint, sector);
1710 }
1b923461 1711 // Double_t value=distPoint[2]-gz;
7f4cb119 1712 if (itype==0){
1713 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1714 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1715 Double_t drphi= r0*(phi1-phi0);
1716 Double_t dr = r1-r0;
1717 distortions[icorr] = drphi-ky*dr;
1718 distortionsR[icorr] = dr;
1719 }
1720 (*pcstream)<<"fit"<<
1721 Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
1722 Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
1723 }
1724 (*pcstream)<<"fit"<<"\n";
1725 }
1726 }
1727 delete pcstream;
1728}
1729
1730
be67055b 1731
b1f0a2a5 1732void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
8b63d99c 1733 //
1734 // make a distortion map out ou fthe residual histogram
1735 // Results are written to the debug streamer - pcstream
1736 // Parameters:
1737 // his0 - input (4D) residual histogram
1738 // pcstream - file to write the tree
1739 // run - run number
1740 // marian.ivanov@cern.ch
1741 const Int_t kMinEntries=50;
1742 Int_t nbins1=his0->GetAxis(1)->GetNbins();
1743 Int_t first1=his0->GetAxis(1)->GetFirst();
1744 Int_t last1 =his0->GetAxis(1)->GetLast();
1745 //
1746 Double_t bz=AliTrackerBase::GetBz();
1747 Int_t idim[4]={0,1,2,3};
1748 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
1749 //
1750 his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1751 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1752 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
1753 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
1754 Int_t first3 = his1->GetAxis(3)->GetFirst();
1755 Int_t last3 = his1->GetAxis(3)->GetLast();
1756 //
1757
1758 for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){ // axis 3 - local angle
1759 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1760 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1761 if (ibin3<first3) {
1762 his1->GetAxis(3)->SetRangeUser(-1,1);
1763 x3=0;
1764 }
1765 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
1766 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1767 Int_t first2 = his3->GetAxis(2)->GetFirst();
1768 Int_t last2 = his3->GetAxis(2)->GetLast();
1769 //
1770 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1771 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1772 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1773 TH1 * hisDelta = his3->Projection(0);
1774 //
1775 Double_t entries = hisDelta->GetEntries();
1776 Double_t mean=0, rms=0;
1777 if (entries>kMinEntries){
1778 mean = hisDelta->GetMean();
1779 rms = hisDelta->GetRMS();
1780 }
1781 (*pcstream)<<hname<<
1782 "run="<<run<<
1783 "bz="<<bz<<
1784 "theta="<<x1<<
1785 "phi="<<x2<<
1786 "snp="<<x3<<
1787 "entries="<<entries<<
1788 "mean="<<mean<<
1789 "rms="<<rms<<
1790 "\n";
1791 delete hisDelta;
1792 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1793 }
1794 delete his3;
1795 }
1796 delete his1;
1797 }
1798}
1799
1800
1801
1802
1803
ffab0c37 1804void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1805 //
1806 // Store object in the OCDB
1807 // By default the object is stored in the current directory
1808 // default comment consit of user name and the date
1809 //
1810 TString ocdbStorage="";
1811 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1812 AliCDBMetaData *metaData= new AliCDBMetaData();
1813 metaData->SetObjectClassName("AliTPCCorrection");
1814 metaData->SetResponsible("Marian Ivanov");
1815 metaData->SetBeamPeriod(1);
1816 metaData->SetAliRootVersion("05-25-01"); //root version
1817 TString userName=gSystem->GetFromPipe("echo $USER");
1818 TString date=gSystem->GetFromPipe("date");
1819
1820 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1821 if (comment) metaData->SetComment(comment);
1822 AliCDBId* id1=NULL;
1823 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1824 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1825 gStorage->Put(this, (*id1), metaData);
1826}
1827
ca58ed4e 1828
7d85e147 1829void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
c9cbd2f2 1830 //
1831 // Fast method to simulate the influence of the given distortion on the vertex reconstruction
1832 //
ca58ed4e 1833
c9cbd2f2 1834 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
1835 if (!magF) AliError("Magneticd field - not initialized");
1836 Double_t bz = magF->SolenoidField(); //field in kGauss
1837 printf("bz: %lf\n",bz);
1838 AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
ca58ed4e 1839
c9cbd2f2 1840 TObjArray aTrk; // Original Track array of Aside
1841 TObjArray daTrk; // Distorted Track array of A side
1842 UShort_t *aId = new UShort_t[nTracks]; // A side Track ID
1843 TObjArray cTrk;
1844 TObjArray dcTrk;
1845 UShort_t *cId = new UShort_t [nTracks];
1846 Int_t id=0;
ca58ed4e 1847 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
7d85e147 1848 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
ca58ed4e 1849 fpt.SetParameters(7.24,0.120);
1850 fpt.SetNpx(10000);
1851 for(Int_t nt=0; nt<nTracks; nt++){
1852 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
7d85e147 1853 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
c9cbd2f2 1854 Double_t pt = fpt.GetRandom(); // momentum for f1
1855 // printf("phi %lf eta %lf pt %lf\n",phi,eta,pt);
ca58ed4e 1856 Short_t sign=1;
1857 if(gRandom->Rndm() < 0.5){
1858 sign =1;
1859 }else{
1860 sign=-1;
1861 }
1862
1863 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1864 Double_t pxyz[3];
1865 pxyz[0]=pt*TMath::Cos(phi);
1866 pxyz[1]=pt*TMath::Sin(phi);
1867 pxyz[2]=pt*TMath::Tan(theta);
1868 Double_t cv[21]={0};
1869 AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
1870
1871 Double_t refX=1.;
1872 Int_t dir=-1;
1873 AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir, NULL);
1874 if (!td) continue;
1875 if (pcstream) (*pcstream)<<"track"<<
1876 "eta="<<eta<<
1877 "theta="<<theta<<
1878 "tOrig.="<<t<<
1879 "td.="<<td<<
1880 "\n";
7d85e147 1881 if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
ca58ed4e 1882 if (td){
c9cbd2f2 1883 daTrk.AddLast(td);
1884 aTrk.AddLast(t);
1885 Int_t nn=aTrk.GetEntriesFast();
1886 aId[nn]=id;
ca58ed4e 1887 }
7d85e147 1888 }else if(( eta<-0.07 )&&( eta>-etaCuts )){
ca58ed4e 1889 if (td){
c9cbd2f2 1890 dcTrk.AddLast(td);
1891 cTrk.AddLast(t);
1892 Int_t nn=cTrk.GetEntriesFast();
1893 cId[nn]=id;
ca58ed4e 1894 }
1895 }
c9cbd2f2 1896 id++;
ca58ed4e 1897 }// end of track loop
1898
1899 vertexer->SetTPCMode();
1900 vertexer->SetConstraintOff();
1901
c9cbd2f2 1902 aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));
1903 avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
1904 cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));
1905 cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
ca58ed4e 1906 if (pcstream) (*pcstream)<<"vertex"<<
1907 "x="<<orgVertex[0]<<
1908 "y="<<orgVertex[1]<<
1909 "z="<<orgVertex[2]<<
1910 "av.="<<&aV<< // distorted vertex A side
1911 "cv.="<<&cV<< // distroted vertex C side
1912 "avO.="<<&avOrg<< // original vertex A side
1913 "cvO.="<<&cvOrg<<
1914 "\n";
c9cbd2f2 1915 delete []aId;
1916 delete []cId;
ca58ed4e 1917}
f1817479 1918
1919void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
1920 //
1921 // make correction available for visualization using
1922 // TFormula, TFX and TTree::Draw
1923 // important in order to check corrections and also compute dervied variables
1924 // e.g correction partial derivatives
1925 //
1926 // NOTE - class is not owner of correction
1927 //
1928 if (!fgVisualCorrection) fgVisualCorrection=new TObjArray;
d0c1b341 1929 if (position!=0&&position>=fgVisualCorrection->GetEntriesFast())
1930 fgVisualCorrection->Expand(position*2);
f1817479 1931 fgVisualCorrection->AddAt(corr, position);
1932}
1933
1934
1935
287fbdfa 1936Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t r, Double_t kZ, Int_t axisType, Int_t corrType){
f1817479 1937 //
1938 // calculate the correction at given position - check the geffCorr
1939 //
1940 if (!fgVisualCorrection) return 0;
1941 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1942 if (!corr) return 0;
25732bff 1943
f1817479 1944 Double_t phi=sector*TMath::Pi()/9.;
287fbdfa 1945 Double_t gx = r*TMath::Cos(phi);
1946 Double_t gy = r*TMath::Sin(phi);
1947 Double_t gz = r*kZ;
f1817479 1948 Int_t nsector=(gz>0) ? 0:18;
1949 //
1950 //
1951 //
1952 Float_t distPoint[3]={gx,gy,gz};
1953 corr->DistortPoint(distPoint, nsector);
1954 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1955 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1956 Double_t phi0=TMath::ATan2(gy,gx);
1957 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1958 if (axisType==0) return r1-r0;
1959 if (axisType==1) return (phi1-phi0)*r0;
1960 if (axisType==2) return distPoint[2]-gz;
1961 return phi1-phi0;
1962}
1963
1964Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
1965 //
1966 // return correction at given x,y,z
1967 //
1968 if (!fgVisualCorrection) return 0;
1969 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1970 if (!corr) return 0;
1971 Double_t phi0= TMath::ATan2(gy,gx);
1972 Int_t nsector=(gz>0) ? 0:18;
1973 Float_t distPoint[3]={gx,gy,gz};
1974 corr->DistortPoint(distPoint, nsector);
1975 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1976 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1977 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1978 if (axisType==0) return r1-r0;
1979 if (axisType==1) return (phi1-phi0)*r0;
1980 if (axisType==2) return distPoint[2]-gz;
1981 return phi1-phi0;
1982}