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