bugfix: adding libAliHLTTPC.so and libAliHLTITS.so which have not been loaded automat...
[u/mrichter/AliRoot.git] / TPC / AliTPCCorrection.cxx
CommitLineData
0116859c 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16////////////////////////////////////////////////////////////////////////////////
17// //
18// AliTPCCorrection class //
19// //
20// This class provides a general framework to deal with space point //
21// distortions. An correction class which inherits from here is for example //
22// AliTPCExBBShape or AliTPCExBTwist //
23// //
24// General functions are (for example): //
25// CorrectPoint(x,roc) where x is the vector of inital positions in //
26// cartesian coordinates and roc represents the Read Out chamber number //
27// according to the offline naming convention. The vector x is overwritten //
28// with the corrected coordinates. //
29// //
30// An alternative usage would be CorrectPoint(x,roc,dx), which leaves the //
31// vector x untouched, put returns the distortions via the vector dx //
32// //
33// The class allows "effective Omega Tau" corrections to be shifted to the //
34// single distortion classes. //
35// //
36// Note: This class is normally used via the class AliTPCComposedCorrection //
37// //
38// date: 27/04/2010 //
39// Authors: Magnus Mager, Stefan Rossegger, Jim Thomas //
40////////////////////////////////////////////////////////////////////////////////
be67055b 41#include "Riostream.h"
0116859c 42
43#include <TH2F.h>
44#include <TMath.h>
45#include <TROOT.h>
cf5b0aa0 46#include <TTreeStream.h>
ffab0c37 47#include <TTree.h>
48#include <TFile.h>
e527a1b9 49#include <TTimeStamp.h>
ffab0c37 50#include <AliCDBStorage.h>
51#include <AliCDBId.h>
52#include <AliCDBMetaData.h>
7f4cb119 53#include "TVectorD.h"
54
1b923461 55
56#include "TRandom.h"
57#include "AliExternalTrackParam.h"
58#include "AliTrackPointArray.h"
59#include "TDatabasePDG.h"
60#include "AliTrackerBase.h"
61#include "AliTPCROC.h"
62#include "THnSparse.h"
63
7f4cb119 64#include "TRandom.h"
65#include "AliTPCTransform.h"
66#include "AliTPCcalibDB.h"
67#include "AliTPCExB.h"
68#include "AliTPCCorrection.h"
69#include "AliTPCRecoParam.h"
cf5b0aa0 70
71#include "AliExternalTrackParam.h"
72#include "AliTrackPointArray.h"
73#include "TDatabasePDG.h"
74#include "AliTrackerBase.h"
75#include "AliTPCROC.h"
8b63d99c 76#include "THnSparse.h"
1b923461 77
7f4cb119 78#include "AliTPCLaserTrack.h"
ca58ed4e 79#include "AliESDVertex.h"
80#include "AliVertexerTracks.h"
81#include "TDatabasePDG.h"
82#include "TF1.h"
0116859c 83
84#include "AliTPCCorrection.h"
1b923461 85#include "AliLog.h"
0116859c 86
cf5b0aa0 87ClassImp(AliTPCCorrection)
88
0116859c 89// FIXME: the following values should come from the database
b1f0a2a5 90const Double_t AliTPCCorrection::fgkTPCZ0 =249.7; // nominal gating grid position
0116859c 91const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
92const Double_t AliTPCCorrection::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
93const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
94const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
95const Double_t AliTPCCorrection::fgkGG =-70.0; // Gating Grid voltage (volts)
96
97
98// FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
99const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
10084.0, 84.5, 85.0, 85.5, 86.0, 87.0, 88.0,
10190.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
102110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
103130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
104150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
105170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
106190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
107210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
108230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 242.0, 244.0, 246.0, 248.0,
109249.0, 249.5, 250.0, 251.5, 252.0 } ;
110
111const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = {
112-249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
113-240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
114-220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
115-200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
116-180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
117-160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
118-140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
119-120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
120-100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0,
121-80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0,
122-60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0,
123-40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0,
124-20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0,
125-1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0,
126 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
127 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0,
128 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0,
129 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0,
130 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0,
131102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
132122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
133142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
134162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
135182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
136202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0,
137222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0,
138242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ;
139
140
141
142AliTPCCorrection::AliTPCCorrection()
534fd34a 143 : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 144{
145 //
146 // default constructor
147 //
148}
149
150AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
534fd34a 151: TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 152{
153 //
154 // default constructor, that set the name and title
155 //
156}
157
158AliTPCCorrection::~AliTPCCorrection() {
159 //
160 // virtual destructor
161 //
162}
163
164void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
165 //
166 // Corrects the initial coordinates x (cartesian coordinates)
167 // according to the given effect (inherited classes)
168 // roc represents the TPC read out chamber (offline numbering convention)
169 //
170 Float_t dx[3];
171 GetCorrection(x,roc,dx);
172 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
173}
174
175void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
176 //
177 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
178 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
179 // roc represents the TPC read out chamber (offline numbering convention)
180 //
181 Float_t dx[3];
182 GetCorrection(x,roc,dx);
183 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
184}
185
186void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
187 //
188 // Distorts the initial coordinates x (cartesian coordinates)
189 // according to the given effect (inherited classes)
190 // roc represents the TPC read out chamber (offline numbering convention)
191 //
192 Float_t dx[3];
193 GetDistortion(x,roc,dx);
194 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
195}
196
197void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
198 //
199 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
200 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
201 // roc represents the TPC read out chamber (offline numbering convention)
202 //
203 Float_t dx[3];
204 GetDistortion(x,roc,dx);
205 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
206}
207
208void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
209 //
210 // This function delivers the correction values dx in respect to the inital coordinates x
211 // roc represents the TPC read out chamber (offline numbering convention)
212 // Note: The dx is overwritten by the inherited effectice class ...
213 //
214 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
215}
216
217void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
218 //
219 // This function delivers the distortion values dx in respect to the inital coordinates x
220 // roc represents the TPC read out chamber (offline numbering convention)
221 //
222 GetCorrection(x,roc,dx);
223 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
224}
225
226void AliTPCCorrection::Init() {
227 //
228 // Initialization funtion (not used at the moment)
229 //
230}
231
e527a1b9 232void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
233 //
234 // Update function
235 //
236}
237
0116859c 238void AliTPCCorrection::Print(Option_t* /*option*/) const {
239 //
240 // Print function to check which correction classes are used
241 // option=="d" prints details regarding the setted magnitude
242 // option=="a" prints the C0 and C1 coefficents for calibration purposes
243 //
244 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
245}
246
534fd34a 247void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
0116859c 248 //
249 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
250 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
251 // calibration run
252 //
534fd34a 253 fT1=t1;
254 fT2=t2;
255 //SetOmegaTauT1T2(omegaTau, t1, t2);
0116859c 256}
257
258TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
259 //
260 // Simple plot functionality.
261 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
262 // in respect to position z within the XY plane.
263 // The histogramm has nx times ny entries.
264 //
265
266 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
267 nx,-250.,250.,ny,-250.,250.);
268 Float_t x[3],dx[3];
269 x[2]=z;
270 Int_t roc=z>0.?0:18; // FIXME
271 for (Int_t iy=1;iy<=ny;++iy) {
272 x[1]=h->GetYaxis()->GetBinCenter(iy);
273 for (Int_t ix=1;ix<=nx;++ix) {
274 x[0]=h->GetXaxis()->GetBinCenter(ix);
275 GetCorrection(x,roc,dx);
276 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
277 if (90.<=r0 && r0<=250.) {
278 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
279 h->SetBinContent(ix,iy,r1-r0);
280 }
281 else
282 h->SetBinContent(ix,iy,0.);
283 }
284 }
285 return h;
286}
287
288TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
289 //
290 // Simple plot functionality.
291 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
292 // in respect to position z within the XY plane.
293 // The histogramm has nx times ny entries.
294 //
295
296 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
297 nx,-250.,250.,ny,-250.,250.);
298 Float_t x[3],dx[3];
299 x[2]=z;
300 Int_t roc=z>0.?0:18; // FIXME
301 for (Int_t iy=1;iy<=ny;++iy) {
302 x[1]=h->GetYaxis()->GetBinCenter(iy);
303 for (Int_t ix=1;ix<=nx;++ix) {
304 x[0]=h->GetXaxis()->GetBinCenter(ix);
305 GetCorrection(x,roc,dx);
306 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
307 if (90.<=r0 && r0<=250.) {
308 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
309 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
310
311 Float_t dphi=phi1-phi0;
312 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
313 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
314
315 h->SetBinContent(ix,iy,r0*dphi);
316 }
317 else
318 h->SetBinContent(ix,iy,0.);
319 }
320 }
321 return h;
322}
323
324TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
325 //
326 // Simple plot functionality.
327 // Returns a 2d hisogram which represents the corrections in r direction (dr)
328 // in respect to angle phi within the ZR plane.
329 // The histogramm has nx times ny entries.
330 //
331 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
332 nz,-250.,250.,nr,85.,250.);
333 Float_t x[3],dx[3];
334 for (Int_t ir=1;ir<=nr;++ir) {
335 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
336 x[0]=radius*TMath::Cos(phi);
337 x[1]=radius*TMath::Sin(phi);
338 for (Int_t iz=1;iz<=nz;++iz) {
339 x[2]=h->GetXaxis()->GetBinCenter(iz);
340 Int_t roc=x[2]>0.?0:18; // FIXME
341 GetCorrection(x,roc,dx);
342 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
343 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
344 h->SetBinContent(iz,ir,r1-r0);
345 }
346 }
0116859c 347 return h;
348
349}
350
351TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
352 //
353 // Simple plot functionality.
354 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
355 // in respect to angle phi within the ZR plane.
356 // The histogramm has nx times ny entries.
357 //
358 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
359 nz,-250.,250.,nr,85.,250.);
360 Float_t x[3],dx[3];
361 for (Int_t iz=1;iz<=nz;++iz) {
362 x[2]=h->GetXaxis()->GetBinCenter(iz);
363 Int_t roc=x[2]>0.?0:18; // FIXME
364 for (Int_t ir=1;ir<=nr;++ir) {
365 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
366 x[0]=radius*TMath::Cos(phi);
367 x[1]=radius*TMath::Sin(phi);
368 GetCorrection(x,roc,dx);
369 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
370 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
371 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
372
373 Float_t dphi=phi1-phi0;
374 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
375 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
376
377 h->SetBinContent(iz,ir,r0*dphi);
378 }
379 }
380 return h;
381}
382
383TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
384 const char *xlabel,const char *ylabel,const char *zlabel,
385 Int_t nbinsx,Double_t xlow,Double_t xup,
386 Int_t nbinsy,Double_t ylow,Double_t yup) {
387 //
388 // Helper function to create a 2d histogramm of given size
389 //
390
391 TString hname=name;
392 Int_t i=0;
393 if (gDirectory) {
394 while (gDirectory->FindObject(hname.Data())) {
395 hname =name;
396 hname+="_";
397 hname+=i;
398 ++i;
399 }
400 }
401 TH2F *h=new TH2F(hname.Data(),title,
402 nbinsx,xlow,xup,
403 nbinsy,ylow,yup);
404 h->GetXaxis()->SetTitle(xlabel);
405 h->GetYaxis()->SetTitle(ylabel);
406 h->GetZaxis()->SetTitle(zlabel);
407 h->SetStats(0);
408 return h;
409}
410
411
412// Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
413
414void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
b1f0a2a5 415 const Double_t er[kNZ][kNR], Double_t &erValue ) {
0116859c 416 //
417 // Interpolate table - 2D interpolation
418 //
b1f0a2a5 419 Double_t saveEr[10] ;
0116859c 420
421 Search( kNZ, fgkZList, z, fJLow ) ;
422 Search( kNR, fgkRList, r, fKLow ) ;
423 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
424 if ( fKLow < 0 ) fKLow = 0 ;
425 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
426 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
427
428 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
b1f0a2a5 429 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
0116859c 430 }
b1f0a2a5 431 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
0116859c 432
433}
434
435
436Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
b1f0a2a5 437 const Int_t order, const Double_t x ) {
0116859c 438 //
439 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
440 //
441
442 Double_t y ;
443 if ( order == 2 ) { // Quadratic Interpolation = 2
444 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
445 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
446 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
447 } else { // Linear Interpolation = 1
448 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
449 }
450
451 return (y);
452
453}
454
455
b1f0a2a5 456void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
0116859c 457 //
458 // Search an ordered table by starting at the most recently used point
459 //
460
461 Long_t middle, high ;
462 Int_t ascend = 0, increment = 1 ;
463
464 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
465
466 if ( low < 0 || low > n-1 ) {
467 low = -1 ; high = n ;
468 } else { // Ordered Search phase
469 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
470 if ( low == n-1 ) return ;
471 high = low + 1 ;
472 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
473 low = high ;
474 increment *= 2 ;
475 high = low + increment ;
476 if ( high > n-1 ) { high = n ; break ; }
477 }
478 } else {
479 if ( low == 0 ) { low = -1 ; return ; }
480 high = low - 1 ;
481 while ( (Int_t)( x < xArray[low] ) == ascend ) {
482 high = low ;
483 increment *= 2 ;
484 if ( increment >= high ) { low = -1 ; break ; }
485 else low = high - increment ;
486 }
487 }
488 }
489
490 while ( (high-low) != 1 ) { // Binary Search Phase
491 middle = ( high + low ) / 2 ;
492 if ( (Int_t)( x >= xArray[middle] ) == ascend )
493 low = middle ;
494 else
495 high = middle ;
496 }
497
498 if ( x == xArray[n-1] ) low = n-2 ;
499 if ( x == xArray[0] ) low = 0 ;
500
501}
502
1b923461 503void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, const TMatrixD &chargeDensity,
504 TMatrixD &arrayErOverEz, const Int_t rows,
505 const Int_t columns, const Int_t iterations ) {
506 //
507 // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
508 //
509 // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
510 // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
511 // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
512 // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
513 // you wish to solve Laplaces equation however it should not contain random numbers or you will get
514 // random numbers back as a solution.
515 // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
516 // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
517 // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
518 // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
519 // number of points and solves the problem again. This happens several times until the maximum number
520 // of points has been included in the array.
521 //
522 // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
523 // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
524 //
525 // Original code by Jim Thomas (STAR TPC Collaboration)
526 //
527
528 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
529
530 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
531 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
532 const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
533
534 TMatrixD arrayEr(rows,columns) ;
535 TMatrixD arrayEz(rows,columns) ;
536
537 //Check that number of rows and columns is suitable for a binary expansion
538
539 if ( !IsPowerOfTwo(rows-1) ) {
540 AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
541 return;
542 }
543 if ( !IsPowerOfTwo(columns-1) ) {
544 AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
545 return;
546 }
547
548 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
549 // Allow for different size grid spacing in R and Z directions
550 // Use a binary expansion of the size of the matrix to speed up the solution of the problem
551
552 Int_t iOne = (rows-1)/4 ;
553 Int_t jOne = (columns-1)/4 ;
554 // Solve for N in 2**N, add one.
555 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
556
557 for ( Int_t count = 0 ; count < loops ; count++ ) {
558 // Loop while the matrix expands & the resolution increases.
559
560 Float_t tempGridSizeR = gridSizeR * iOne ;
561 Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
562 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
563
564 // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
565 std::vector<float> coef1(rows) ;
566 std::vector<float> coef2(rows) ;
567
568 for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
569 Float_t radius = fgkIFCRadius + i*gridSizeR ;
570 coef1[i] = 1.0 + tempGridSizeR/(2*radius);
571 coef2[i] = 1.0 - tempGridSizeR/(2*radius);
572 }
573
574 TMatrixD sumChargeDensity(rows,columns) ;
575
576 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
577 Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
578 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
579 if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
580 else {
581 // Add up all enclosed charge density contributions within 1/2 unit in all directions
582 Float_t weight = 0.0 ;
583 Float_t sum = 0.0 ;
584 sumChargeDensity(i,j) = 0.0 ;
585 for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
586 for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
587 if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
588 else
589 weight = 1.0 ;
590 // Note that this is cylindrical geometry
591 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
592 sum += weight*radius ;
593 }
594 }
595 sumChargeDensity(i,j) /= sum ;
596 }
597 sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
598 }
599 }
600
601 for ( Int_t k = 1 ; k <= iterations; k++ ) {
602 // Solve Poisson's Equation
603 // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
604 // as interations increase.
605 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
606 Float_t overRelaxM1 = overRelax - 1.0 ;
607 Float_t overRelaxtempFourth, overRelaxcoef5 ;
608 overRelaxtempFourth = overRelax * tempFourth ;
609 overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
610
611 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
612 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
613
614 arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
615 + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
616 - overRelaxcoef5 * arrayV(i,j)
617 + coef1[i] * arrayV(i+iOne,j)
618 + sumChargeDensity(i,j)
619 ) * overRelaxtempFourth;
620 }
621 }
622
623 if ( k == iterations ) {
624 // After full solution is achieved, copy low resolution solution into higher res array
625 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
626 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
627
628 if ( iOne > 1 ) {
629 arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
630 if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
631 }
632 if ( jOne > 1 ) {
633 arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
634 if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
635 }
636 if ( iOne > 1 && jOne > 1 ) {
637 arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
638 if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
639 if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
640 // Note that this leaves a point at the upper left and lower right corners uninitialized.
641 // -> Not a big deal.
642 }
643
644 }
645 }
646 }
647
648 }
649
650 iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
651 jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
652
653 }
654
655 // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
656 for ( Int_t j = 0 ; j < columns ; j++ ) {
657 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
658 arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
659 arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
660 }
661
662 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
663 for ( Int_t i = 0 ; i < rows ; i++) {
664 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
665 arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
666 arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
667 }
668
669 for ( Int_t i = 0 ; i < rows ; i++) {
670 // Note: go back and compare to old version of this code. See notes below.
671 // JT Test ... attempt to divide by real Ez not Ez to first order
672 for ( Int_t j = 0 ; j < columns ; j++ ) {
673 arrayEz(i,j) += ezField;
674 // This adds back the overall Z gradient of the field (main E field component)
675 }
676 // Warning: (-=) assumes you are using an error potetial without the overall Field included
677 }
678
679 // Integrate Er/Ez from Z to zero
680 for ( Int_t j = 0 ; j < columns ; j++ ) {
681 for ( Int_t i = 0 ; i < rows ; i++ ) {
682 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
683 arrayErOverEz(i,j) = 0.0 ;
684 for ( Int_t k = j ; k < columns ; k++ ) {
685 arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
686 if ( index != 4 ) index = 4; else index = 2 ;
687 }
688 if ( index == 4 ) arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
689 if ( index == 2 ) arrayErOverEz(i,j) +=
690 (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
691 -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
692 if ( j == columns-2 ) arrayErOverEz(i,j) =
693 (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
694 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
695 if ( j == columns-1 ) arrayErOverEz(i,j) = 0.0 ;
696 }
697 }
698
699}
700
701
702
710bda39 703Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1b923461 704 //
705 // Helperfunction: Check if integer is a power of 2
706 //
707 Int_t j = 0;
708 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
709 if ( j == 1 ) return(1) ; // True
710 return(0) ; // False
711}
712
cf5b0aa0 713
b1f0a2a5 714AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
cf5b0aa0 715 //
716 // Fit the track parameters - without and with distortion
717 // 1. Space points in the TPC are simulated along the trajectory
718 // 2. Space points distorted
719 // 3. Fits the non distorted and distroted track to the reference plane at refX
720 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
721 //
722 // trackIn - input track parameters
723 // refX - reference X to fit the track
724 // dir - direction - out=1 or in=-1
725 // pcstream - debug streamer to check the results
726 //
cad404e1 727 // see AliExternalTrackParam.h documentation:
728 // track1.fP[0] - local y (rphi)
729 // track1.fP[1] - z
730 // track1.fP[2] - sinus of local inclination angle
731 // track1.fP[3] - tangent of deep angle
732 // track1.fP[4] - 1/pt
1b923461 733
cf5b0aa0 734 AliTPCROC * roc = AliTPCROC::Instance();
735 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
736 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
737 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
cf5b0aa0 738 const Double_t kMaxSnp = 0.85;
739 const Double_t kSigmaY=0.1;
740 const Double_t kSigmaZ=0.1;
ca58ed4e 741 const Double_t kMaxR=500;
742 const Double_t kMaxZ=500;
cf5b0aa0 743 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
ca58ed4e 744 Int_t npoints1=0;
745 Int_t npoints2=0;
cf5b0aa0 746
be67055b 747 AliExternalTrackParam track(trackIn); //
cf5b0aa0 748 // generate points
749 AliTrackPointArray pointArray0(npoints0);
750 AliTrackPointArray pointArray1(npoints0);
751 Double_t xyz[3];
ca58ed4e 752 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 753 //
754 // simulate the track
755 Int_t npoints=0;
756 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
757 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
ca58ed4e 758 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 759 track.GetXYZ(xyz);
ca58ed4e 760 xyz[0]+=gRandom->Gaus(0,0.00005);
761 xyz[1]+=gRandom->Gaus(0,0.00005);
762 xyz[2]+=gRandom->Gaus(0,0.00005);
763 if (TMath::Abs(track.GetZ())>kMaxZ) break;
764 if (TMath::Abs(track.GetX())>kMaxR) break;
cf5b0aa0 765 AliTrackPoint pIn0; // space point
766 AliTrackPoint pIn1;
ffab0c37 767 Int_t sector= (xyz[2]>0)? 0:18;
cf5b0aa0 768 pointArray0.GetPoint(pIn0,npoints);
769 pointArray1.GetPoint(pIn1,npoints);
770 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
771 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
ffab0c37 772 DistortPoint(distPoint, sector);
cf5b0aa0 773 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
774 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
775 //
776 track.Rotate(alpha);
777 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
778 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
779 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
780 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
781 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
782 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
783 pointArray0.AddPoint(npoints, &pIn0);
784 pointArray1.AddPoint(npoints, &pIn1);
785 npoints++;
786 if (npoints>=npoints0) break;
787 }
7f4cb119 788 if (npoints<npoints0/2) return 0;
cf5b0aa0 789 //
790 // refit track
791 //
792 AliExternalTrackParam *track0=0;
793 AliExternalTrackParam *track1=0;
794 AliTrackPoint point1,point2,point3;
795 if (dir==1) { //make seed inner
796 pointArray0.GetPoint(point1,1);
4486a91f 797 pointArray0.GetPoint(point2,30);
798 pointArray0.GetPoint(point3,60);
cf5b0aa0 799 }
800 if (dir==-1){ //make seed outer
4486a91f 801 pointArray0.GetPoint(point1,npoints-60);
802 pointArray0.GetPoint(point2,npoints-30);
cf5b0aa0 803 pointArray0.GetPoint(point3,npoints-1);
804 }
805 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
806 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
807
cf5b0aa0 808 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
8b63d99c 809 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
cf5b0aa0 810 //
811 AliTrackPoint pIn0;
812 AliTrackPoint pIn1;
813 pointArray0.GetPoint(pIn0,ipoint);
814 pointArray1.GetPoint(pIn1,ipoint);
815 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
816 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
817 //
ca58ed4e 818 if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
819 if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
820 if (TMath::Abs(track0->GetZ())>kMaxZ) break;
821 if (TMath::Abs(track0->GetX())>kMaxR) break;
822 if (TMath::Abs(track1->GetZ())>kMaxZ) break;
823 if (TMath::Abs(track1->GetX())>kMaxR) break;
824
8b63d99c 825 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
cf5b0aa0 826 //
827 Double_t pointPos[2]={0,0};
828 Double_t pointCov[3]={0,0,0};
829 pointPos[0]=prot0.GetY();//local y
830 pointPos[1]=prot0.GetZ();//local z
831 pointCov[0]=prot0.GetCov()[3];//simay^2
832 pointCov[1]=prot0.GetCov()[4];//sigmayz
833 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
ca58ed4e 834 if (!track0->Update(pointPos,pointCov)) break;
cf5b0aa0 835 //
8b63d99c 836 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
837 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
838 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
839
0b736a46 840 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
841 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
cf5b0aa0 842 pointCov[0]=prot1.GetCov()[3];//simay^2
843 pointCov[1]=prot1.GetCov()[4];//sigmayz
844 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
ca58ed4e 845 if (!track1->Update(pointPos,pointCov)) break;
846 npoints1++;
847 npoints2++;
cf5b0aa0 848 }
ca58ed4e 849 if (npoints2<npoints) return 0;
8b63d99c 850 AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
cf5b0aa0 851 track1->Rotate(track0->GetAlpha());
ca58ed4e 852 AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
cf5b0aa0 853
cad404e1 854 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
cf5b0aa0 855 "point0.="<<&pointArray0<< // points
856 "point1.="<<&pointArray1<< // distorted points
857 "trackIn.="<<&track<< // original track
858 "track0.="<<track0<< // fitted track
859 "track1.="<<track1<< // fitted distorted track
860 "\n";
be67055b 861 new(&trackIn) AliExternalTrackParam(*track0);
cf5b0aa0 862 delete track0;
863 return track1;
864}
865
866
ffab0c37 867
868
869
870TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
871 //
872 // create the distortion tree on a mesh with granularity given by step
873 // return the tree with distortions at given position
874 // Map is created on the mesh with given step size
875 //
876 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
877 Float_t xyz[3];
878 for (Double_t x= -250; x<250; x+=step){
879 for (Double_t y= -250; y<250; y+=step){
880 Double_t r = TMath::Sqrt(x*x+y*y);
881 if (r<80) continue;
882 if (r>250) continue;
883 for (Double_t z= -250; z<250; z+=step){
884 Int_t roc=(z>0)?0:18;
885 xyz[0]=x;
886 xyz[1]=y;
887 xyz[2]=z;
888 Double_t phi = TMath::ATan2(y,x);
889 DistortPoint(xyz,roc);
890 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
891 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
892 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
893 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
894 Double_t dx = xyz[0]-x;
895 Double_t dy = xyz[1]-y;
896 Double_t dz = xyz[2]-z;
897 Double_t dr=r1-r;
898 Double_t drphi=(phi1-phi)*r;
899 (*pcstream)<<"distortion"<<
900 "x="<<x<< // original position
901 "y="<<y<<
902 "z="<<z<<
903 "r="<<r<<
904 "phi="<<phi<<
905 "x1="<<xyz[0]<< // distorted position
906 "y1="<<xyz[1]<<
907 "z1="<<xyz[2]<<
908 "r1="<<r1<<
909 "phi1="<<phi1<<
910 //
911 "dx="<<dx<< // delta position
912 "dy="<<dy<<
913 "dz="<<dz<<
914 "dr="<<dr<<
915 "drphi="<<drphi<<
916 "\n";
917 }
918 }
919 }
920 delete pcstream;
921 TFile f(Form("correction%s.root",GetName()));
922 TTree * tree = (TTree*)f.Get("distortion");
923 TTree * tree2= tree->CopyTree("1");
924 tree2->SetName(Form("dist%s",GetName()));
925 tree2->SetDirectory(0);
926 delete tree;
927 return tree2;
928}
929
930
931
be67055b 932
b1f0a2a5 933void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
be67055b 934 //
935 // Make a fit tree:
936 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
937 // calculates partial distortions
938 // Partial distortion is stored in the resulting tree
939 // Output is storred in the file distortion_<dettype>_<partype>.root
940 // Partial distortion is stored with the name given by correction name
941 //
942 //
943 // Parameters of function:
944 // input - input tree
945 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex
946 // ppype - parameter type
947 // corrArray - array with partial corrections
948 // step - skipe entries - if 1 all entries processed - it is slow
949 // debug 0 if debug on also space points dumped - it is slow
b322e06a 950 const Double_t kMaxSnp = 0.85;
951 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
952 // const Double_t kB2C=-0.299792458e-3;
7f4cb119 953 const Int_t kMinEntries=50;
be67055b 954 Double_t phi,theta, snp, mean,rms, entries;
955 tinput->SetBranchAddress("theta",&theta);
956 tinput->SetBranchAddress("phi", &phi);
957 tinput->SetBranchAddress("snp",&snp);
958 tinput->SetBranchAddress("mean",&mean);
959 tinput->SetBranchAddress("rms",&rms);
960 tinput->SetBranchAddress("entries",&entries);
961 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
962 //
963 Int_t nentries=tinput->GetEntries();
964 Int_t ncorr=corrArray->GetEntries();
7f4cb119 965 Double_t corrections[100]={0}; //
be67055b 966 Double_t tPar[5];
967 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
968 Double_t refX=0;
969 Int_t dir=0;
7f4cb119 970 if (dtype==0) {refX=85.; dir=-1;}
971 if (dtype==1) {refX=275.; dir=1;}
972 if (dtype==2) {refX=85.; dir=-1;}
973 if (dtype==3) {refX=360.; dir=-1;}
be67055b 974 //
975 for (Int_t ientry=0; ientry<nentries; ientry+=step){
976 tinput->GetEntry(ientry);
7f4cb119 977 if (TMath::Abs(snp)>kMaxSnp) continue;
be67055b 978 tPar[0]=0;
979 tPar[1]=theta*refX;
980 tPar[2]=snp;
981 tPar[3]=theta;
4486a91f 982 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
8b63d99c 983 Double_t bz=AliTrackerBase::GetBz();
4486a91f 984 if (refX>10. && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
985 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
7f4cb119 986 AliExternalTrackParam track(refX,phi,tPar,cov);
987 Double_t xyz[3];
988 track.GetXYZ(xyz);
989 Int_t id=0;
990 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
0b736a46 991 if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature
be67055b 992 (*pcstream)<<"fit"<<
8b63d99c 993 "bz="<<bz<< // magnetic filed used
be67055b 994 "dtype="<<dtype<< // detector match type
995 "ptype="<<ptype<< // parameter type
996 "theta="<<theta<< // theta
997 "phi="<<phi<< // phi
998 "snp="<<snp<< // snp
999 "mean="<<mean<< // mean dist value
1000 "rms="<<rms<< // rms
7f4cb119 1001 "gx="<<xyz[0]<< // global position at reference
1002 "gy="<<xyz[1]<< // global position at reference
1003 "gz="<<xyz[2]<< // global position at reference
1004 "dRrec="<<dRrec<< // delta Radius in reconstruction
1005 "id="<<id<< // track id
be67055b 1006 "entries="<<entries;// number of entries in bin
1007 //
1008 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1009 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1010 corrections[icorr]=0;
1011 if (entries>kMinEntries){
1012 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1013 AliExternalTrackParam *trackOut = 0;
1014 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1015 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
7f4cb119 1016 if (dtype==0) {refX=85.; dir=-1;}
1017 if (dtype==1) {refX=275.; dir=1;}
b1f0a2a5 1018 if (dtype==2) {refX=0; dir=-1;}
7f4cb119 1019 if (dtype==3) {refX=360.; dir=-1;}
b1f0a2a5 1020 //
7f4cb119 1021 if (trackOut){
1022 AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1023 trackOut->Rotate(trackIn.GetAlpha());
1024 trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1025 //
1026 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1027 delete trackOut;
1028 }else{
1029 corrections[icorr]=0;
1030 }
0b736a46 1031 if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature
be67055b 1032 }
7f4cb119 1033 Double_t dRdummy=0;
be67055b 1034 (*pcstream)<<"fit"<<
7f4cb119 1035 Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
1036 Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
1037 // for points it is neccessary
be67055b 1038 }
1039 (*pcstream)<<"fit"<<"\n";
1040 }
1041 delete pcstream;
1042}
1043
1044
1045
7f4cb119 1046void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1047 //
1048 // Make a laser fit tree for global minimization
1049 //
1050 const Double_t cutErrY=0.1;
1051 const Double_t cutErrZ=0.1;
1052 const Double_t kEpsilon=0.00000001;
1053 TVectorD *vecdY=0;
1054 TVectorD *vecdZ=0;
1055 TVectorD *veceY=0;
1056 TVectorD *veceZ=0;
1057 AliTPCLaserTrack *ltr=0;
1058 AliTPCLaserTrack::LoadTracks();
1059 tree->SetBranchAddress("dY.",&vecdY);
1060 tree->SetBranchAddress("dZ.",&vecdZ);
1061 tree->SetBranchAddress("eY.",&veceY);
1062 tree->SetBranchAddress("eZ.",&veceZ);
1063 tree->SetBranchAddress("LTr.",&ltr);
1064 Int_t entries= tree->GetEntries();
1065 TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1066 Double_t bz=AliTrackerBase::GetBz();
1067 //
1068
1069 for (Int_t ientry=0; ientry<entries; ientry++){
1070 tree->GetEntry(ientry);
1071 if (!ltr->GetVecGX()){
1072 ltr->UpdatePoints();
1073 }
1074 TVectorD * delta= (itype==0)? vecdY:vecdZ;
1075 TVectorD * err= (itype==0)? veceY:veceZ;
1076
1077 for (Int_t irow=0; irow<159; irow++){
1078 Int_t nentries = 1000;
1079 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1080 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1081 Int_t dtype=4;
1082 Double_t phi =(*ltr->GetVecPhi())[irow];
1083 Double_t theta =ltr->GetTgl();
1084 Double_t mean=delta->GetMatrixArray()[irow];
1085 Double_t gx=0,gy=0,gz=0;
1086 Double_t snp = (*ltr->GetVecP2())[irow];
1087 Double_t rms = 0.1+err->GetMatrixArray()[irow];
1088 gx = (*ltr->GetVecGX())[irow];
1089 gy = (*ltr->GetVecGY())[irow];
1090 gz = (*ltr->GetVecGZ())[irow];
1091 Int_t bundle= ltr->GetBundle();
1092 Double_t dRrec=0;
1093 //
1094 // get delta R used in reconstruction
1095 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
1096 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1097 const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1098 Double_t xyz0[3]={gx,gy,gz};
1099 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1100 //
1101 // old ExB correction
1102 //
1103 if(recoParam&&recoParam->GetUseExBCorrection()) {
1104 Double_t xyz1[3]={gx,gy,gz};
1105 calib->GetExB()->Correct(xyz0,xyz1);
1106 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1107 dRrec=oldR-newR;
1108 }
1109 if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1110 Float_t xyz1[3]={gx,gy,gz};
1111 Int_t sector=(gz>0)?0:18;
1112 correction->CorrectPoint(xyz1, sector);
1113 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1114 dRrec=oldR-newR;
1115 }
1116
1117
1118 (*pcstream)<<"fit"<<
1119 "bz="<<bz<< // magnetic filed used
1120 "dtype="<<dtype<< // detector match type
1121 "ptype="<<itype<< // parameter type
1122 "theta="<<theta<< // theta
1123 "phi="<<phi<< // phi
1124 "snp="<<snp<< // snp
1125 "mean="<<mean<< // mean dist value
1126 "rms="<<rms<< // rms
1127 "gx="<<gx<< // global position
1128 "gy="<<gy<< // global position
1129 "gz="<<gz<< // global position
1130 "dRrec="<<dRrec<< // delta Radius in reconstruction
1131 "id="<<bundle<< //bundle
1132 "entries="<<nentries;// number of entries in bin
1133 //
1134 //
1135 Double_t ky = TMath::Tan(TMath::ASin(snp));
1136 Int_t ncorr = corrArray->GetEntries();
1137 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
1138 Double_t phi0 = TMath::ATan2(gy,gx);
1139 Double_t distortions[1000]={0};
1140 Double_t distortionsR[1000]={0};
1141 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1142 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1143 Float_t distPoint[3]={gx,gy,gz};
1144 Int_t sector= (gz>0)? 0:18;
1145 if (r0>80){
1146 corr->DistortPoint(distPoint, sector);
1147 }
1b923461 1148 // Double_t value=distPoint[2]-gz;
7f4cb119 1149 if (itype==0){
1150 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1151 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1152 Double_t drphi= r0*(phi1-phi0);
1153 Double_t dr = r1-r0;
1154 distortions[icorr] = drphi-ky*dr;
1155 distortionsR[icorr] = dr;
1156 }
1157 (*pcstream)<<"fit"<<
1158 Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
1159 Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
1160 }
1161 (*pcstream)<<"fit"<<"\n";
1162 }
1163 }
1164 delete pcstream;
1165}
1166
1167
be67055b 1168
b1f0a2a5 1169void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
8b63d99c 1170 //
1171 // make a distortion map out ou fthe residual histogram
1172 // Results are written to the debug streamer - pcstream
1173 // Parameters:
1174 // his0 - input (4D) residual histogram
1175 // pcstream - file to write the tree
1176 // run - run number
1177 // marian.ivanov@cern.ch
1178 const Int_t kMinEntries=50;
1179 Int_t nbins1=his0->GetAxis(1)->GetNbins();
1180 Int_t first1=his0->GetAxis(1)->GetFirst();
1181 Int_t last1 =his0->GetAxis(1)->GetLast();
1182 //
1183 Double_t bz=AliTrackerBase::GetBz();
1184 Int_t idim[4]={0,1,2,3};
1185 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
1186 //
1187 his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1188 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1189 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
1190 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
1191 Int_t first3 = his1->GetAxis(3)->GetFirst();
1192 Int_t last3 = his1->GetAxis(3)->GetLast();
1193 //
1194
1195 for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){ // axis 3 - local angle
1196 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1197 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1198 if (ibin3<first3) {
1199 his1->GetAxis(3)->SetRangeUser(-1,1);
1200 x3=0;
1201 }
1202 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
1203 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1204 Int_t first2 = his3->GetAxis(2)->GetFirst();
1205 Int_t last2 = his3->GetAxis(2)->GetLast();
1206 //
1207 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1208 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1209 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1210 TH1 * hisDelta = his3->Projection(0);
1211 //
1212 Double_t entries = hisDelta->GetEntries();
1213 Double_t mean=0, rms=0;
1214 if (entries>kMinEntries){
1215 mean = hisDelta->GetMean();
1216 rms = hisDelta->GetRMS();
1217 }
1218 (*pcstream)<<hname<<
1219 "run="<<run<<
1220 "bz="<<bz<<
1221 "theta="<<x1<<
1222 "phi="<<x2<<
1223 "snp="<<x3<<
1224 "entries="<<entries<<
1225 "mean="<<mean<<
1226 "rms="<<rms<<
1227 "\n";
1228 delete hisDelta;
1229 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1230 }
1231 delete his3;
1232 }
1233 delete his1;
1234 }
1235}
1236
1237
1238
1239
1240
ffab0c37 1241void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1242 //
1243 // Store object in the OCDB
1244 // By default the object is stored in the current directory
1245 // default comment consit of user name and the date
1246 //
1247 TString ocdbStorage="";
1248 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1249 AliCDBMetaData *metaData= new AliCDBMetaData();
1250 metaData->SetObjectClassName("AliTPCCorrection");
1251 metaData->SetResponsible("Marian Ivanov");
1252 metaData->SetBeamPeriod(1);
1253 metaData->SetAliRootVersion("05-25-01"); //root version
1254 TString userName=gSystem->GetFromPipe("echo $USER");
1255 TString date=gSystem->GetFromPipe("date");
1256
1257 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1258 if (comment) metaData->SetComment(comment);
1259 AliCDBId* id1=NULL;
1260 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1261 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1262 gStorage->Put(this, (*id1), metaData);
1263}
1264
ca58ed4e 1265
7d85e147 1266void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
ca58ed4e 1267
1268 AliVertexerTracks *vertexer = new AliVertexerTracks(5);// 5kGaus
1269
1270 TObjArray ATrk; // Original Track array of Aside
1271 TObjArray dATrk; // Distorted Track array of A side
1272 UShort_t *AId = new UShort_t[nTracks]; // A side Track ID
1273 TObjArray CTrk;
1274 TObjArray dCTrk;
1275 UShort_t *CId = new UShort_t [nTracks];
1276 Int_t ID=0;
1277 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
7d85e147 1278 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
ca58ed4e 1279 fpt.SetParameters(7.24,0.120);
1280 fpt.SetNpx(10000);
1281 for(Int_t nt=0; nt<nTracks; nt++){
1282 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
7d85e147 1283 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
ca58ed4e 1284 Double_t pt = fpt.GetRandom();// momentum for f1
1285 Short_t sign=1;
1286 if(gRandom->Rndm() < 0.5){
1287 sign =1;
1288 }else{
1289 sign=-1;
1290 }
1291
1292 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1293 Double_t pxyz[3];
1294 pxyz[0]=pt*TMath::Cos(phi);
1295 pxyz[1]=pt*TMath::Sin(phi);
1296 pxyz[2]=pt*TMath::Tan(theta);
1297 Double_t cv[21]={0};
1298 AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
1299
1300 Double_t refX=1.;
1301 Int_t dir=-1;
1302 AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir, NULL);
1303 if (!td) continue;
1304 if (pcstream) (*pcstream)<<"track"<<
1305 "eta="<<eta<<
1306 "theta="<<theta<<
1307 "tOrig.="<<t<<
1308 "td.="<<td<<
1309 "\n";
7d85e147 1310 if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
ca58ed4e 1311 if (td){
1312 dATrk.AddLast(td);
1313 ATrk.AddLast(t);
1314 Int_t nn=ATrk.GetEntriesFast();
1315 AId[nn]=ID;
1316 }
7d85e147 1317 }else if(( eta<-0.07 )&&( eta>-etaCuts )){
ca58ed4e 1318 if (td){
1319 dCTrk.AddLast(td);
1320 CTrk.AddLast(t);
1321 Int_t nn=CTrk.GetEntriesFast();
1322 CId[nn]=ID;
1323 }
1324 }
1325 ID++;
1326 }// end of track loop
1327
1328 vertexer->SetTPCMode();
1329 vertexer->SetConstraintOff();
1330
1331 aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dATrk,AId));
1332 avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&ATrk,AId));
1333 cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dCTrk,CId));
1334 cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&CTrk,CId));
1335 if (pcstream) (*pcstream)<<"vertex"<<
1336 "x="<<orgVertex[0]<<
1337 "y="<<orgVertex[1]<<
1338 "z="<<orgVertex[2]<<
1339 "av.="<<&aV<< // distorted vertex A side
1340 "cv.="<<&cV<< // distroted vertex C side
1341 "avO.="<<&avOrg<< // original vertex A side
1342 "cvO.="<<&cvOrg<<
1343 "\n";
1344 delete []AId;
1345 delete []CId;
1346}