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