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