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