]> 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
fixing some minor coverty reports
[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 //
b1f0a2a5 503 Double_t saveEr[10] ;
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
526 Double_t saveEr[10], savedEr[10] ;
527 Double_t saveEphi[10], savedEphi[10] ;
528 Double_t saveEz[10], savedEz[10] ;
529
530 Search( kNZ, fgkZList, z, fILow ) ;
531 Search( kNPhi, fgkPhiList, z, fJLow ) ;
532 Search( kNR, fgkRList, r, fKLow ) ;
533
534 if ( fILow < 0 ) fILow = 0 ; // check if out of range
535 if ( fJLow < 0 ) fJLow = 0 ;
536 if ( fKLow < 0 ) fKLow = 0 ;
537
538 if ( fILow + order >= kNZ - 1 ) fILow = kNZ - 1 - order ;
539 if ( fJLow + order >= kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
540 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
541
542 for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
543 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
544 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r ) ;
545 saveEphi[j-fJLow] = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
546 saveEz[j-fJLow] = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r ) ;
547 }
548 savedEr[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi ) ;
549 savedEphi[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
550 savedEz[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi ) ;
551 }
552 erValue = Interpolate( &fgkZList[fILow], savedEr, order, z ) ;
553 ephiValue = Interpolate( &fgkZList[fILow], savedEphi, order, z ) ;
554 ezValue = Interpolate( &fgkZList[fILow], savedEz, order, z ) ;
555
556}
557
558Double_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
559 const Int_t nx, const Int_t ny, const Double_t xv[], const Double_t yv[],
560 const TMatrixD &array ) {
561 //
562 // Interpolate table (TMatrix format) - 2D interpolation
563 //
564
565 static Int_t jlow = 0, klow = 0 ;
566 Double_t saveArray[10] ;
567
568 Search( nx, xv, x, jlow ) ;
569 Search( ny, yv, y, klow ) ;
570 if ( jlow < 0 ) jlow = 0 ; // check if out of range
571 if ( klow < 0 ) klow = 0 ;
572 if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
573 if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
574
575 for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
576 {
577 Double_t *ajkl = &((TMatrixD&)array)(j,klow);
578 saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
579 }
580
581 return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
582
583}
584
585Double_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x, const Double_t y, const Double_t z,
586 const Int_t nx, const Int_t ny, const Int_t nz,
587 const Double_t xv[], const Double_t yv[], const Double_t zv[],
588 TMatrixD **arrayofArrays ) {
589 //
590 // Interpolate table (TMatrix format) - 3D interpolation
591 //
592
593 static Int_t ilow = 0, jlow = 0, klow = 0 ;
594 Double_t saveArray[10], savedArray[10] ;
595
596 Search( nx, xv, x, ilow ) ;
597 Search( ny, yv, y, jlow ) ;
598 Search( nz, zv, z, klow ) ;
599
600 if ( ilow < 0 ) ilow = 0 ; // check if out of range
601 if ( jlow < 0 ) jlow = 0 ;
602 if ( klow < 0 ) klow = 0 ;
603
604 if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
605 if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
606 if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
607
608 for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
609 {
610 TMatrixD &table = *arrayofArrays[k] ;
611 for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
612 {
613 saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
614 }
615 savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
616 }
617 return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
618
619}
620
621
0116859c 622
623Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
b1f0a2a5 624 const Int_t order, const Double_t x ) {
0116859c 625 //
626 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
627 //
628
629 Double_t y ;
630 if ( order == 2 ) { // Quadratic Interpolation = 2
631 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
632 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
633 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
634 } else { // Linear Interpolation = 1
635 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
636 }
637
638 return (y);
639
640}
641
642
b1f0a2a5 643void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
0116859c 644 //
645 // Search an ordered table by starting at the most recently used point
646 //
647
648 Long_t middle, high ;
649 Int_t ascend = 0, increment = 1 ;
650
651 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
652
653 if ( low < 0 || low > n-1 ) {
654 low = -1 ; high = n ;
655 } else { // Ordered Search phase
656 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
657 if ( low == n-1 ) return ;
658 high = low + 1 ;
659 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
660 low = high ;
661 increment *= 2 ;
662 high = low + increment ;
663 if ( high > n-1 ) { high = n ; break ; }
664 }
665 } else {
666 if ( low == 0 ) { low = -1 ; return ; }
667 high = low - 1 ;
668 while ( (Int_t)( x < xArray[low] ) == ascend ) {
669 high = low ;
670 increment *= 2 ;
671 if ( increment >= high ) { low = -1 ; break ; }
672 else low = high - increment ;
673 }
674 }
675 }
676
677 while ( (high-low) != 1 ) { // Binary Search Phase
678 middle = ( high + low ) / 2 ;
679 if ( (Int_t)( x >= xArray[middle] ) == ascend )
680 low = middle ;
681 else
682 high = middle ;
683 }
684
685 if ( x == xArray[n-1] ) low = n-2 ;
686 if ( x == xArray[0] ) low = 0 ;
687
688}
689
c9cbd2f2 690void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
691 TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
692 const Int_t rows, const Int_t columns, const Int_t iterations,
693 const Bool_t rocDisplacement ) {
1b923461 694 //
695 // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
696 //
697 // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
698 // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
699 // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
700 // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
701 // you wish to solve Laplaces equation however it should not contain random numbers or you will get
702 // random numbers back as a solution.
703 // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
704 // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
705 // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
706 // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
707 // number of points and solves the problem again. This happens several times until the maximum number
708 // of points has been included in the array.
709 //
710 // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
711 // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
712 //
c9cbd2f2 713 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
714 //
1b923461 715 // Original code by Jim Thomas (STAR TPC Collaboration)
716 //
717
718 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
719
720 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
721 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
722 const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
723
724 TMatrixD arrayEr(rows,columns) ;
725 TMatrixD arrayEz(rows,columns) ;
726
727 //Check that number of rows and columns is suitable for a binary expansion
728
729 if ( !IsPowerOfTwo(rows-1) ) {
730 AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
731 return;
732 }
733 if ( !IsPowerOfTwo(columns-1) ) {
734 AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
735 return;
736 }
737
738 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
739 // Allow for different size grid spacing in R and Z directions
740 // Use a binary expansion of the size of the matrix to speed up the solution of the problem
741
742 Int_t iOne = (rows-1)/4 ;
743 Int_t jOne = (columns-1)/4 ;
744 // Solve for N in 2**N, add one.
745 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
746
747 for ( Int_t count = 0 ; count < loops ; count++ ) {
748 // Loop while the matrix expands & the resolution increases.
749
750 Float_t tempGridSizeR = gridSizeR * iOne ;
751 Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
752 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
753
754 // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
755 std::vector<float> coef1(rows) ;
756 std::vector<float> coef2(rows) ;
757
758 for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
759 Float_t radius = fgkIFCRadius + i*gridSizeR ;
760 coef1[i] = 1.0 + tempGridSizeR/(2*radius);
761 coef2[i] = 1.0 - tempGridSizeR/(2*radius);
762 }
763
764 TMatrixD sumChargeDensity(rows,columns) ;
765
766 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
767 Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
768 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
769 if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
770 else {
771 // Add up all enclosed charge density contributions within 1/2 unit in all directions
772 Float_t weight = 0.0 ;
773 Float_t sum = 0.0 ;
774 sumChargeDensity(i,j) = 0.0 ;
775 for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
776 for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
777 if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
778 else
779 weight = 1.0 ;
780 // Note that this is cylindrical geometry
781 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
782 sum += weight*radius ;
783 }
784 }
785 sumChargeDensity(i,j) /= sum ;
786 }
787 sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
788 }
789 }
790
791 for ( Int_t k = 1 ; k <= iterations; k++ ) {
792 // Solve Poisson's Equation
793 // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
794 // as interations increase.
795 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
796 Float_t overRelaxM1 = overRelax - 1.0 ;
797 Float_t overRelaxtempFourth, overRelaxcoef5 ;
798 overRelaxtempFourth = overRelax * tempFourth ;
799 overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
800
801 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
802 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
803
804 arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
805 + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
806 - overRelaxcoef5 * arrayV(i,j)
807 + coef1[i] * arrayV(i+iOne,j)
808 + sumChargeDensity(i,j)
809 ) * overRelaxtempFourth;
810 }
811 }
812
813 if ( k == iterations ) {
814 // After full solution is achieved, copy low resolution solution into higher res array
815 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
816 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
817
818 if ( iOne > 1 ) {
819 arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
820 if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
821 }
822 if ( jOne > 1 ) {
823 arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
824 if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
825 }
826 if ( iOne > 1 && jOne > 1 ) {
827 arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
828 if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
829 if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
830 // Note that this leaves a point at the upper left and lower right corners uninitialized.
831 // -> Not a big deal.
832 }
833
834 }
835 }
836 }
837
838 }
839
840 iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
841 jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
842
c9cbd2f2 843 sumChargeDensity.Clear();
1b923461 844 }
845
846 // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
847 for ( Int_t j = 0 ; j < columns ; j++ ) {
848 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
849 arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
850 arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
851 }
852
853 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
854 for ( Int_t i = 0 ; i < rows ; i++) {
855 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
856 arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
857 arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
858 }
859
860 for ( Int_t i = 0 ; i < rows ; i++) {
861 // Note: go back and compare to old version of this code. See notes below.
862 // JT Test ... attempt to divide by real Ez not Ez to first order
863 for ( Int_t j = 0 ; j < columns ; j++ ) {
864 arrayEz(i,j) += ezField;
865 // This adds back the overall Z gradient of the field (main E field component)
866 }
867 // Warning: (-=) assumes you are using an error potetial without the overall Field included
868 }
869
870 // Integrate Er/Ez from Z to zero
871 for ( Int_t j = 0 ; j < columns ; j++ ) {
872 for ( Int_t i = 0 ; i < rows ; i++ ) {
c9cbd2f2 873
1b923461 874 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
875 arrayErOverEz(i,j) = 0.0 ;
c9cbd2f2 876 arrayDeltaEz(i,j) = 0.0 ;
877
1b923461 878 for ( Int_t k = j ; k < columns ; k++ ) {
879 arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
c9cbd2f2 880 arrayDeltaEz(i,j) += index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
1b923461 881 if ( index != 4 ) index = 4; else index = 2 ;
882 }
c9cbd2f2 883 if ( index == 4 ) {
884 arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
885 arrayDeltaEz(i,j) -= (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
886 }
887 if ( index == 2 ) {
888 arrayErOverEz(i,j) += (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
889 -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
890 arrayDeltaEz(i,j) += (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
891 -2.5*(arrayEz(i,columns-1)-ezField));
892 }
893 if ( j == columns-2 ) {
894 arrayErOverEz(i,j) = (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
895 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
896 arrayDeltaEz(i,j) = (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
897 +1.5*(arrayEz(i,columns-1)-ezField) ) ;
898 }
899 if ( j == columns-1 ) {
900 arrayErOverEz(i,j) = 0.0 ;
901 arrayDeltaEz(i,j) = 0.0 ;
902 }
1b923461 903 }
904 }
905
c9cbd2f2 906 // calculate z distortion from the integrated Delta Ez residuals
907 // and include the aquivalence (Volt to cm) of the ROC shift !!
908
909 for ( Int_t j = 0 ; j < columns ; j++ ) {
910 for ( Int_t i = 0 ; i < rows ; i++ ) {
911
912 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
913 arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
914
915 // ROC Potential in cm aquivalent
916 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
917 if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift; // add the ROC misaligment
918
919 }
920 }
921
922 arrayEr.Clear();
923 arrayEz.Clear();
924
1b923461 925}
926
c9cbd2f2 927void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
928 TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
929 const Int_t rows, const Int_t columns, const Int_t phislices,
930 const Float_t deltaphi, const Int_t iterations, const Int_t symmetry,
931 Bool_t rocDisplacement ) {
932 //
933 // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
934 //
935 // NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
936 // The number of rows and COLUMNS can be different.
937 //
938 // ROWS == 2**M + 1
939 // COLUMNS == 2**N + 1
940 // PHISLICES == Arbitrary but greater than 3
941 //
942 // DeltaPhi in Radians
943 //
944 // SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
945 // = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
946 //
947 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
948
949 const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
950
951 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
952 const Float_t gridSizePhi = deltaphi ;
953 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
954 const Float_t ratioPhi = gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
955 const Float_t ratioZ = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
956
957 TMatrixD arrayE(rows,columns) ;
958
959 // Check that the number of rows and columns is suitable for a binary expansion
960 if ( !IsPowerOfTwo((rows-1)) ) {
961 AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
962 return; }
963 if ( !IsPowerOfTwo((columns-1)) ) {
964 AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
965 return; }
966 if ( phislices <= 3 ) {
967 AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
968 return; }
969 if ( phislices > 1000 ) {
970 AliError("Poisson3D phislices > 1000 is not allowed (nor wise) ");
971 return; }
972
973 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
974 // Allow for different size grid spacing in R and Z directions
975 // Use a binary expansion of the matrix to speed up the solution of the problem
976
977 Int_t loops, mplus, mminus, signplus, signminus ;
978 Int_t ione = (rows-1)/4 ;
979 Int_t jone = (columns-1)/4 ;
980 loops = TMath::Max(ione, jone) ; // Calculate the number of loops for the binary expansion
981 loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ; // Solve for N in 2**N
982
983 TMatrixD* arrayofSumChargeDensities[1000] ; // Create temporary arrays to store low resolution charge arrays
984
985 for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
986
987 for ( Int_t count = 0 ; count < loops ; count++ ) { // START the master loop and do the binary expansion
988
989 Float_t tempgridSizeR = gridSizeR * ione ;
990 Float_t tempratioPhi = ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
991 Float_t tempratioZ = ratioZ * ione * ione / ( jone * jone ) ;
992
993 std::vector<float> coef1(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
994 std::vector<float> coef2(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
995 std::vector<float> coef3(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
996 std::vector<float> coef4(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
997
998 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
999 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1000 coef1[i] = 1.0 + tempgridSizeR/(2*radius);
1001 coef2[i] = 1.0 - tempgridSizeR/(2*radius);
1002 coef3[i] = tempratioPhi/(radius*radius);
1003 coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
1004 }
1005
1006 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1007 TMatrixD &chargeDensity = *arrayofChargeDensities[m] ;
1008 TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
1009 for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
1010 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1011 for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
1012 if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1013 else { // Add up all enclosed charge density contributions within 1/2 unit in all directions
1014 Float_t weight = 0.0 ;
1015 Float_t sum = 0.0 ;
1016 sumChargeDensity(i,j) = 0.0 ;
1017 for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
1018 for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
1019 if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
1020 else
1021 weight = 1.0 ;
1022 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1023 sum += weight*radius ;
1024 }
1025 }
1026 sumChargeDensity(i,j) /= sum ;
1027 }
1028 sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
1029 }
1030 }
1031 }
1032
1033 for ( Int_t k = 1 ; k <= iterations; k++ ) {
1034
1035 // over-relaxation index, >= 1 but < 2
1036 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1037 Float_t overRelaxM1 = overRelax - 1.0 ;
1038
1039 std::vector<float> overRelaxcoef4(rows) ; // Do this the standard C++ way to avoid gcc extensions
1040 std::vector<float> overRelaxcoef5(rows) ; // Do this the standard C++ way to avoid gcc extensions
1041
1042 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1043 overRelaxcoef4[i] = overRelax * coef4[i] ;
1044 overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
1045 }
1046
1047 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1048
1049 mplus = m + 1; signplus = 1 ;
1050 mminus = m - 1 ; signminus = 1 ;
1051 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1052 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1053 if ( mminus < 0 ) mminus = 1 ;
1054 }
1055 else if (symmetry==-1) { // Anti-symmetry in phi
1056 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1057 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1058 }
1059 else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1060 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1061 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1062 }
1063 TMatrixD& arrayV = *arrayofArrayV[m] ;
1064 TMatrixD& arrayVP = *arrayofArrayV[mplus] ;
1065 TMatrixD& arrayVM = *arrayofArrayV[mminus] ;
1066 TMatrixD& sumChargeDensity = *arrayofSumChargeDensities[m] ;
1067
1068 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1069 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1070
1071 arrayV(i,j) = ( coef2[i] * arrayV(i-ione,j)
1072 + tempratioZ * ( arrayV(i,j-jone) + arrayV(i,j+jone) )
1073 - overRelaxcoef5[i] * arrayV(i,j)
1074 + coef1[i] * arrayV(i+ione,j)
1075 + coef3[i] * ( signplus*arrayVP(i,j) + signminus*arrayVM(i,j) )
1076 + sumChargeDensity(i,j)
1077 ) * overRelaxcoef4[i] ;
1078 // Note: over-relax the solution at each step. This speeds up the convergance.
1079
1080 }
1081 }
1082
1083 if ( k == iterations ) { // After full solution is achieved, copy low resolution solution into higher res array
1084 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1085 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1086
1087 if ( ione > 1 ) {
1088 arrayV(i+ione/2,j) = ( arrayV(i+ione,j) + arrayV(i,j) ) / 2 ;
1089 if ( i == ione ) arrayV(i-ione/2,j) = ( arrayV(0,j) + arrayV(ione,j) ) / 2 ;
1090 }
1091 if ( jone > 1 ) {
1092 arrayV(i,j+jone/2) = ( arrayV(i,j+jone) + arrayV(i,j) ) / 2 ;
1093 if ( j == jone ) arrayV(i,j-jone/2) = ( arrayV(i,0) + arrayV(i,jone) ) / 2 ;
1094 }
1095 if ( ione > 1 && jone > 1 ) {
1096 arrayV(i+ione/2,j+jone/2) = ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
1097 if ( i == ione ) arrayV(i-ione/2,j-jone/2) = ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
1098 if ( j == jone ) arrayV(i-ione/2,j-jone/2) = ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
1099 // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
1100 }
1101 }
1102 }
1103 }
1104
1105 }
1106 }
1107
1108 ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
1109 jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
1110
1111 }
1112
1113 //Differentiate V(r) and solve for E(r) using special equations for the first and last row
1114 //Integrate E(r)/E(z) from point of origin to pad plane
1115
1116 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1117 TMatrixD& arrayV = *arrayofArrayV[m] ;
1118 TMatrixD& eroverEz = *arrayofEroverEz[m] ;
1119
1120 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1121
1122 // Differentiate in R
1123 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1124 arrayE(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1125 arrayE(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1126 // Integrate over Z
1127 for ( Int_t i = 0 ; i < rows ; i++ ) {
1128 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1129 eroverEz(i,j) = 0.0 ;
1130 for ( Int_t k = j ; k < columns ; k++ ) {
1131
1132 eroverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1133 if ( index != 4 ) index = 4; else index = 2 ;
1134 }
1135 if ( index == 4 ) eroverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1136 if ( index == 2 ) eroverEz(i,j) +=
1137 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1138 if ( j == columns-2 ) eroverEz(i,j) =
1139 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1140 if ( j == columns-1 ) eroverEz(i,j) = 0.0 ;
1141 }
1142 }
1143 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1144 // eroverEz.Draw("surf") ; } // JT test
1145 }
1146
1147 //Differentiate V(r) and solve for E(phi)
1148 //Integrate E(phi)/E(z) from point of origin to pad plane
1149
1150 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1151
1152 mplus = m + 1; signplus = 1 ;
1153 mminus = m - 1 ; signminus = 1 ;
1154 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1155 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1156 if ( mminus < 0 ) mminus = 1 ;
1157 }
1158 else if (symmetry==-1) { // Anti-symmetry in phi
1159 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1160 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1161 }
1162 else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
1163 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1164 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1165 }
1166 TMatrixD &arrayVP = *arrayofArrayV[mplus] ;
1167 TMatrixD &arrayVM = *arrayofArrayV[mminus] ;
1168 TMatrixD &ePhioverEz = *arrayofEPhioverEz[m] ;
1169 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1170 // Differentiate in Phi
1171 for ( Int_t i = 0 ; i < rows ; i++ ) {
1172 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1173 arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
1174 }
1175 // Integrate over Z
1176 for ( Int_t i = 0 ; i < rows ; i++ ) {
1177 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1178 ePhioverEz(i,j) = 0.0 ;
1179 for ( Int_t k = j ; k < columns ; k++ ) {
1180
1181 ePhioverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1182 if ( index != 4 ) index = 4; else index = 2 ;
1183 }
1184 if ( index == 4 ) ePhioverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1185 if ( index == 2 ) ePhioverEz(i,j) +=
1186 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1187 if ( j == columns-2 ) ePhioverEz(i,j) =
1188 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1189 if ( j == columns-1 ) ePhioverEz(i,j) = 0.0 ;
1190 }
1191 }
1192 // if ( m == 5 ) { TCanvas* c2 = new TCanvas("arrayE","arrayE",50,50,840,600) ; c2 -> cd() ;
1193 // arrayE.Draw("surf") ; } // JT test
1194 }
1195
1196
1197 // Differentiate V(r) and solve for E(z) using special equations for the first and last row
1198 // Integrate (E(z)-Ezstd) from point of origin to pad plane
1199
1200 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1201 TMatrixD& arrayV = *arrayofArrayV[m] ;
1202 TMatrixD& deltaEz = *arrayofDeltaEz[m] ;
1203
1204 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1205 for ( Int_t i = 0 ; i < rows ; i++) {
1206 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1207 arrayE(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1208 arrayE(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1209 }
1210
1211 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1212 // Integrate over Z
1213 for ( Int_t i = 0 ; i < rows ; i++ ) {
1214 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1215 deltaEz(i,j) = 0.0 ;
1216 for ( Int_t k = j ; k < columns ; k++ ) {
1217 deltaEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k) ;
1218 if ( index != 4 ) index = 4; else index = 2 ;
1219 }
1220 if ( index == 4 ) deltaEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1) ;
1221 if ( index == 2 ) deltaEz(i,j) +=
1222 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
1223 if ( j == columns-2 ) deltaEz(i,j) =
1224 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
1225 if ( j == columns-1 ) deltaEz(i,j) = 0.0 ;
1226 }
1227 }
1228 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1229 // eroverEz.Draw("surf") ; } // JT test
1230
1231 // calculate z distortion from the integrated Delta Ez residuals
1232 // and include the aquivalence (Volt to cm) of the ROC shift !!
1233
1234 for ( Int_t j = 0 ; j < columns ; j++ ) {
1235 for ( Int_t i = 0 ; i < rows ; i++ ) {
1236
1237 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1238 deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
1239
1240 // ROC Potential in cm aquivalent
1241 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
1242 if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift; // add the ROC misaligment
1243
1244 }
1245 }
1246
1247 } // end loop over phi
1248
1249
1250
1251 for ( Int_t k = 0 ; k < phislices ; k++ )
1252 {
1253 arrayofSumChargeDensities[k]->Delete() ;
1254 }
1255
1256
1257
1258 arrayE.Clear();
1259}
1b923461 1260
1261
710bda39 1262Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1b923461 1263 //
1264 // Helperfunction: Check if integer is a power of 2
1265 //
1266 Int_t j = 0;
1267 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
1268 if ( j == 1 ) return(1) ; // True
1269 return(0) ; // False
1270}
1271
cf5b0aa0 1272
b1f0a2a5 1273AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
cf5b0aa0 1274 //
1275 // Fit the track parameters - without and with distortion
1276 // 1. Space points in the TPC are simulated along the trajectory
1277 // 2. Space points distorted
1278 // 3. Fits the non distorted and distroted track to the reference plane at refX
1279 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
1280 //
1281 // trackIn - input track parameters
1282 // refX - reference X to fit the track
1283 // dir - direction - out=1 or in=-1
1284 // pcstream - debug streamer to check the results
1285 //
cad404e1 1286 // see AliExternalTrackParam.h documentation:
1287 // track1.fP[0] - local y (rphi)
1288 // track1.fP[1] - z
1289 // track1.fP[2] - sinus of local inclination angle
1290 // track1.fP[3] - tangent of deep angle
1291 // track1.fP[4] - 1/pt
1b923461 1292
cf5b0aa0 1293 AliTPCROC * roc = AliTPCROC::Instance();
1294 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
1295 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
1296 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
cf5b0aa0 1297 const Double_t kMaxSnp = 0.85;
1298 const Double_t kSigmaY=0.1;
1299 const Double_t kSigmaZ=0.1;
ca58ed4e 1300 const Double_t kMaxR=500;
1301 const Double_t kMaxZ=500;
cf5b0aa0 1302 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
ca58ed4e 1303 Int_t npoints1=0;
1304 Int_t npoints2=0;
cf5b0aa0 1305
be67055b 1306 AliExternalTrackParam track(trackIn); //
cf5b0aa0 1307 // generate points
1308 AliTrackPointArray pointArray0(npoints0);
1309 AliTrackPointArray pointArray1(npoints0);
1310 Double_t xyz[3];
ca58ed4e 1311 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 1312 //
1313 // simulate the track
1314 Int_t npoints=0;
1315 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
1316 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
ca58ed4e 1317 if (!AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 1318 track.GetXYZ(xyz);
ca58ed4e 1319 xyz[0]+=gRandom->Gaus(0,0.00005);
1320 xyz[1]+=gRandom->Gaus(0,0.00005);
1321 xyz[2]+=gRandom->Gaus(0,0.00005);
1322 if (TMath::Abs(track.GetZ())>kMaxZ) break;
1323 if (TMath::Abs(track.GetX())>kMaxR) break;
cf5b0aa0 1324 AliTrackPoint pIn0; // space point
1325 AliTrackPoint pIn1;
ffab0c37 1326 Int_t sector= (xyz[2]>0)? 0:18;
cf5b0aa0 1327 pointArray0.GetPoint(pIn0,npoints);
1328 pointArray1.GetPoint(pIn1,npoints);
1329 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
1330 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
ffab0c37 1331 DistortPoint(distPoint, sector);
cf5b0aa0 1332 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
1333 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
1334 //
1335 track.Rotate(alpha);
1336 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
1337 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
1338 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
1339 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
1340 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
1341 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
1342 pointArray0.AddPoint(npoints, &pIn0);
1343 pointArray1.AddPoint(npoints, &pIn1);
1344 npoints++;
1345 if (npoints>=npoints0) break;
1346 }
7f4cb119 1347 if (npoints<npoints0/2) return 0;
cf5b0aa0 1348 //
1349 // refit track
1350 //
1351 AliExternalTrackParam *track0=0;
1352 AliExternalTrackParam *track1=0;
1353 AliTrackPoint point1,point2,point3;
1354 if (dir==1) { //make seed inner
1355 pointArray0.GetPoint(point1,1);
4486a91f 1356 pointArray0.GetPoint(point2,30);
1357 pointArray0.GetPoint(point3,60);
cf5b0aa0 1358 }
1359 if (dir==-1){ //make seed outer
4486a91f 1360 pointArray0.GetPoint(point1,npoints-60);
1361 pointArray0.GetPoint(point2,npoints-30);
cf5b0aa0 1362 pointArray0.GetPoint(point3,npoints-1);
1363 }
1364 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
1365 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
1366
cf5b0aa0 1367 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
8b63d99c 1368 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
cf5b0aa0 1369 //
1370 AliTrackPoint pIn0;
1371 AliTrackPoint pIn1;
1372 pointArray0.GetPoint(pIn0,ipoint);
1373 pointArray1.GetPoint(pIn1,ipoint);
1374 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
1375 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
1376 //
ca58ed4e 1377 if (!AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1378 if (!AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp)) break;
1379 if (TMath::Abs(track0->GetZ())>kMaxZ) break;
1380 if (TMath::Abs(track0->GetX())>kMaxR) break;
1381 if (TMath::Abs(track1->GetZ())>kMaxZ) break;
1382 if (TMath::Abs(track1->GetX())>kMaxR) break;
1383
8b63d99c 1384 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
cf5b0aa0 1385 //
1386 Double_t pointPos[2]={0,0};
1387 Double_t pointCov[3]={0,0,0};
1388 pointPos[0]=prot0.GetY();//local y
1389 pointPos[1]=prot0.GetZ();//local z
1390 pointCov[0]=prot0.GetCov()[3];//simay^2
1391 pointCov[1]=prot0.GetCov()[4];//sigmayz
1392 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
ca58ed4e 1393 if (!track0->Update(pointPos,pointCov)) break;
cf5b0aa0 1394 //
8b63d99c 1395 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
1396 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
1397 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
1398
0b736a46 1399 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
1400 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
cf5b0aa0 1401 pointCov[0]=prot1.GetCov()[3];//simay^2
1402 pointCov[1]=prot1.GetCov()[4];//sigmayz
1403 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
ca58ed4e 1404 if (!track1->Update(pointPos,pointCov)) break;
1405 npoints1++;
1406 npoints2++;
cf5b0aa0 1407 }
ca58ed4e 1408 if (npoints2<npoints) return 0;
8b63d99c 1409 AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
cf5b0aa0 1410 track1->Rotate(track0->GetAlpha());
ca58ed4e 1411 AliTrackerBase::PropagateTrackToBxByBz(track1,refX,kMass,2.,kTRUE,kMaxSnp);
cf5b0aa0 1412
cad404e1 1413 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
cf5b0aa0 1414 "point0.="<<&pointArray0<< // points
1415 "point1.="<<&pointArray1<< // distorted points
1416 "trackIn.="<<&track<< // original track
1417 "track0.="<<track0<< // fitted track
1418 "track1.="<<track1<< // fitted distorted track
1419 "\n";
be67055b 1420 new(&trackIn) AliExternalTrackParam(*track0);
cf5b0aa0 1421 delete track0;
1422 return track1;
1423}
1424
1425
ffab0c37 1426
1427
1428
1429TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
1430 //
1431 // create the distortion tree on a mesh with granularity given by step
1432 // return the tree with distortions at given position
1433 // Map is created on the mesh with given step size
1434 //
1435 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
1436 Float_t xyz[3];
1437 for (Double_t x= -250; x<250; x+=step){
1438 for (Double_t y= -250; y<250; y+=step){
1439 Double_t r = TMath::Sqrt(x*x+y*y);
1440 if (r<80) continue;
1441 if (r>250) continue;
1442 for (Double_t z= -250; z<250; z+=step){
1443 Int_t roc=(z>0)?0:18;
1444 xyz[0]=x;
1445 xyz[1]=y;
1446 xyz[2]=z;
1447 Double_t phi = TMath::ATan2(y,x);
1448 DistortPoint(xyz,roc);
1449 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1450 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
1451 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
1452 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
1453 Double_t dx = xyz[0]-x;
1454 Double_t dy = xyz[1]-y;
1455 Double_t dz = xyz[2]-z;
1456 Double_t dr=r1-r;
1457 Double_t drphi=(phi1-phi)*r;
1458 (*pcstream)<<"distortion"<<
1459 "x="<<x<< // original position
1460 "y="<<y<<
1461 "z="<<z<<
1462 "r="<<r<<
1463 "phi="<<phi<<
1464 "x1="<<xyz[0]<< // distorted position
1465 "y1="<<xyz[1]<<
1466 "z1="<<xyz[2]<<
1467 "r1="<<r1<<
1468 "phi1="<<phi1<<
1469 //
1470 "dx="<<dx<< // delta position
1471 "dy="<<dy<<
1472 "dz="<<dz<<
1473 "dr="<<dr<<
1474 "drphi="<<drphi<<
1475 "\n";
1476 }
1477 }
1478 }
1479 delete pcstream;
1480 TFile f(Form("correction%s.root",GetName()));
1481 TTree * tree = (TTree*)f.Get("distortion");
1482 TTree * tree2= tree->CopyTree("1");
1483 tree2->SetName(Form("dist%s",GetName()));
1484 tree2->SetDirectory(0);
1485 delete tree;
1486 return tree2;
1487}
1488
1489
1490
be67055b 1491
b1f0a2a5 1492void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
be67055b 1493 //
1494 // Make a fit tree:
1495 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1496 // calculates partial distortions
1497 // Partial distortion is stored in the resulting tree
1498 // Output is storred in the file distortion_<dettype>_<partype>.root
1499 // Partial distortion is stored with the name given by correction name
1500 //
1501 //
1502 // Parameters of function:
1503 // input - input tree
1504 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex
1505 // ppype - parameter type
1506 // corrArray - array with partial corrections
1507 // step - skipe entries - if 1 all entries processed - it is slow
1508 // debug 0 if debug on also space points dumped - it is slow
c9cbd2f2 1509
b322e06a 1510 const Double_t kMaxSnp = 0.85;
1511 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1512 // const Double_t kB2C=-0.299792458e-3;
7f4cb119 1513 const Int_t kMinEntries=50;
be67055b 1514 Double_t phi,theta, snp, mean,rms, entries;
1515 tinput->SetBranchAddress("theta",&theta);
1516 tinput->SetBranchAddress("phi", &phi);
1517 tinput->SetBranchAddress("snp",&snp);
1518 tinput->SetBranchAddress("mean",&mean);
1519 tinput->SetBranchAddress("rms",&rms);
1520 tinput->SetBranchAddress("entries",&entries);
1521 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
1522 //
1523 Int_t nentries=tinput->GetEntries();
1524 Int_t ncorr=corrArray->GetEntries();
7f4cb119 1525 Double_t corrections[100]={0}; //
be67055b 1526 Double_t tPar[5];
1527 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1528 Double_t refX=0;
1529 Int_t dir=0;
7f4cb119 1530 if (dtype==0) {refX=85.; dir=-1;}
1531 if (dtype==1) {refX=275.; dir=1;}
1532 if (dtype==2) {refX=85.; dir=-1;}
1533 if (dtype==3) {refX=360.; dir=-1;}
be67055b 1534 //
1535 for (Int_t ientry=0; ientry<nentries; ientry+=step){
1536 tinput->GetEntry(ientry);
7f4cb119 1537 if (TMath::Abs(snp)>kMaxSnp) continue;
be67055b 1538 tPar[0]=0;
1539 tPar[1]=theta*refX;
1540 tPar[2]=snp;
1541 tPar[3]=theta;
4486a91f 1542 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
8b63d99c 1543 Double_t bz=AliTrackerBase::GetBz();
4486a91f 1544 if (refX>10. && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
1545 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
7f4cb119 1546 AliExternalTrackParam track(refX,phi,tPar,cov);
1547 Double_t xyz[3];
1548 track.GetXYZ(xyz);
1549 Int_t id=0;
1550 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
0b736a46 1551 if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature
be67055b 1552 (*pcstream)<<"fit"<<
8b63d99c 1553 "bz="<<bz<< // magnetic filed used
be67055b 1554 "dtype="<<dtype<< // detector match type
1555 "ptype="<<ptype<< // parameter type
1556 "theta="<<theta<< // theta
1557 "phi="<<phi<< // phi
1558 "snp="<<snp<< // snp
1559 "mean="<<mean<< // mean dist value
1560 "rms="<<rms<< // rms
7f4cb119 1561 "gx="<<xyz[0]<< // global position at reference
1562 "gy="<<xyz[1]<< // global position at reference
1563 "gz="<<xyz[2]<< // global position at reference
1564 "dRrec="<<dRrec<< // delta Radius in reconstruction
1565 "id="<<id<< // track id
be67055b 1566 "entries="<<entries;// number of entries in bin
1567 //
1568 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1569 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1570 corrections[icorr]=0;
1571 if (entries>kMinEntries){
1572 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1573 AliExternalTrackParam *trackOut = 0;
1574 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1575 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
7f4cb119 1576 if (dtype==0) {refX=85.; dir=-1;}
1577 if (dtype==1) {refX=275.; dir=1;}
b1f0a2a5 1578 if (dtype==2) {refX=0; dir=-1;}
7f4cb119 1579 if (dtype==3) {refX=360.; dir=-1;}
b1f0a2a5 1580 //
7f4cb119 1581 if (trackOut){
1582 AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
1583 trackOut->Rotate(trackIn.GetAlpha());
1584 trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1585 //
1586 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1587 delete trackOut;
1588 }else{
1589 corrections[icorr]=0;
1590 }
0b736a46 1591 if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature
be67055b 1592 }
7f4cb119 1593 Double_t dRdummy=0;
be67055b 1594 (*pcstream)<<"fit"<<
7f4cb119 1595 Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
1596 Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
1597 // for points it is neccessary
be67055b 1598 }
1599 (*pcstream)<<"fit"<<"\n";
1600 }
1601 delete pcstream;
1602}
1603
1604
1605
7f4cb119 1606void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
1607 //
1608 // Make a laser fit tree for global minimization
1609 //
1610 const Double_t cutErrY=0.1;
1611 const Double_t cutErrZ=0.1;
1612 const Double_t kEpsilon=0.00000001;
1613 TVectorD *vecdY=0;
1614 TVectorD *vecdZ=0;
1615 TVectorD *veceY=0;
1616 TVectorD *veceZ=0;
1617 AliTPCLaserTrack *ltr=0;
1618 AliTPCLaserTrack::LoadTracks();
1619 tree->SetBranchAddress("dY.",&vecdY);
1620 tree->SetBranchAddress("dZ.",&vecdZ);
1621 tree->SetBranchAddress("eY.",&veceY);
1622 tree->SetBranchAddress("eZ.",&veceZ);
1623 tree->SetBranchAddress("LTr.",&ltr);
1624 Int_t entries= tree->GetEntries();
1625 TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
1626 Double_t bz=AliTrackerBase::GetBz();
1627 //
1628
1629 for (Int_t ientry=0; ientry<entries; ientry++){
1630 tree->GetEntry(ientry);
1631 if (!ltr->GetVecGX()){
1632 ltr->UpdatePoints();
1633 }
1634 TVectorD * delta= (itype==0)? vecdY:vecdZ;
1635 TVectorD * err= (itype==0)? veceY:veceZ;
1636
1637 for (Int_t irow=0; irow<159; irow++){
1638 Int_t nentries = 1000;
1639 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
1640 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
1641 Int_t dtype=4;
1642 Double_t phi =(*ltr->GetVecPhi())[irow];
1643 Double_t theta =ltr->GetTgl();
1644 Double_t mean=delta->GetMatrixArray()[irow];
1645 Double_t gx=0,gy=0,gz=0;
1646 Double_t snp = (*ltr->GetVecP2())[irow];
1647 Double_t rms = 0.1+err->GetMatrixArray()[irow];
1648 gx = (*ltr->GetVecGX())[irow];
1649 gy = (*ltr->GetVecGY())[irow];
1650 gz = (*ltr->GetVecGZ())[irow];
1651 Int_t bundle= ltr->GetBundle();
1652 Double_t dRrec=0;
1653 //
1654 // get delta R used in reconstruction
1655 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
1656 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
1657 const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
1658 Double_t xyz0[3]={gx,gy,gz};
1659 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
1660 //
1661 // old ExB correction
1662 //
1663 if(recoParam&&recoParam->GetUseExBCorrection()) {
1664 Double_t xyz1[3]={gx,gy,gz};
1665 calib->GetExB()->Correct(xyz0,xyz1);
1666 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1667 dRrec=oldR-newR;
1668 }
1669 if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
1670 Float_t xyz1[3]={gx,gy,gz};
1671 Int_t sector=(gz>0)?0:18;
1672 correction->CorrectPoint(xyz1, sector);
1673 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
1674 dRrec=oldR-newR;
1675 }
1676
1677
1678 (*pcstream)<<"fit"<<
1679 "bz="<<bz<< // magnetic filed used
1680 "dtype="<<dtype<< // detector match type
1681 "ptype="<<itype<< // parameter type
1682 "theta="<<theta<< // theta
1683 "phi="<<phi<< // phi
1684 "snp="<<snp<< // snp
1685 "mean="<<mean<< // mean dist value
1686 "rms="<<rms<< // rms
1687 "gx="<<gx<< // global position
1688 "gy="<<gy<< // global position
1689 "gz="<<gz<< // global position
1690 "dRrec="<<dRrec<< // delta Radius in reconstruction
1691 "id="<<bundle<< //bundle
1692 "entries="<<nentries;// number of entries in bin
1693 //
1694 //
1695 Double_t ky = TMath::Tan(TMath::ASin(snp));
1696 Int_t ncorr = corrArray->GetEntries();
1697 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
1698 Double_t phi0 = TMath::ATan2(gy,gx);
1699 Double_t distortions[1000]={0};
1700 Double_t distortionsR[1000]={0};
1701 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1702 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1703 Float_t distPoint[3]={gx,gy,gz};
1704 Int_t sector= (gz>0)? 0:18;
1705 if (r0>80){
1706 corr->DistortPoint(distPoint, sector);
1707 }
1b923461 1708 // Double_t value=distPoint[2]-gz;
7f4cb119 1709 if (itype==0){
1710 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1711 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
1712 Double_t drphi= r0*(phi1-phi0);
1713 Double_t dr = r1-r0;
1714 distortions[icorr] = drphi-ky*dr;
1715 distortionsR[icorr] = dr;
1716 }
1717 (*pcstream)<<"fit"<<
1718 Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
1719 Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
1720 }
1721 (*pcstream)<<"fit"<<"\n";
1722 }
1723 }
1724 delete pcstream;
1725}
1726
1727
be67055b 1728
b1f0a2a5 1729void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
8b63d99c 1730 //
1731 // make a distortion map out ou fthe residual histogram
1732 // Results are written to the debug streamer - pcstream
1733 // Parameters:
1734 // his0 - input (4D) residual histogram
1735 // pcstream - file to write the tree
1736 // run - run number
1737 // marian.ivanov@cern.ch
1738 const Int_t kMinEntries=50;
1739 Int_t nbins1=his0->GetAxis(1)->GetNbins();
1740 Int_t first1=his0->GetAxis(1)->GetFirst();
1741 Int_t last1 =his0->GetAxis(1)->GetLast();
1742 //
1743 Double_t bz=AliTrackerBase::GetBz();
1744 Int_t idim[4]={0,1,2,3};
1745 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
1746 //
1747 his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1748 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
1749 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
1750 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
1751 Int_t first3 = his1->GetAxis(3)->GetFirst();
1752 Int_t last3 = his1->GetAxis(3)->GetLast();
1753 //
1754
1755 for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){ // axis 3 - local angle
1756 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
1757 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
1758 if (ibin3<first3) {
1759 his1->GetAxis(3)->SetRangeUser(-1,1);
1760 x3=0;
1761 }
1762 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
1763 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1764 Int_t first2 = his3->GetAxis(2)->GetFirst();
1765 Int_t last2 = his3->GetAxis(2)->GetLast();
1766 //
1767 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
1768 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
1769 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1770 TH1 * hisDelta = his3->Projection(0);
1771 //
1772 Double_t entries = hisDelta->GetEntries();
1773 Double_t mean=0, rms=0;
1774 if (entries>kMinEntries){
1775 mean = hisDelta->GetMean();
1776 rms = hisDelta->GetRMS();
1777 }
1778 (*pcstream)<<hname<<
1779 "run="<<run<<
1780 "bz="<<bz<<
1781 "theta="<<x1<<
1782 "phi="<<x2<<
1783 "snp="<<x3<<
1784 "entries="<<entries<<
1785 "mean="<<mean<<
1786 "rms="<<rms<<
1787 "\n";
1788 delete hisDelta;
1789 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
1790 }
1791 delete his3;
1792 }
1793 delete his1;
1794 }
1795}
1796
1797
1798
1799
1800
ffab0c37 1801void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1802 //
1803 // Store object in the OCDB
1804 // By default the object is stored in the current directory
1805 // default comment consit of user name and the date
1806 //
1807 TString ocdbStorage="";
1808 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1809 AliCDBMetaData *metaData= new AliCDBMetaData();
1810 metaData->SetObjectClassName("AliTPCCorrection");
1811 metaData->SetResponsible("Marian Ivanov");
1812 metaData->SetBeamPeriod(1);
1813 metaData->SetAliRootVersion("05-25-01"); //root version
1814 TString userName=gSystem->GetFromPipe("echo $USER");
1815 TString date=gSystem->GetFromPipe("date");
1816
1817 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1818 if (comment) metaData->SetComment(comment);
1819 AliCDBId* id1=NULL;
1820 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1821 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1822 gStorage->Put(this, (*id1), metaData);
1823}
1824
ca58ed4e 1825
7d85e147 1826void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
c9cbd2f2 1827 //
1828 // Fast method to simulate the influence of the given distortion on the vertex reconstruction
1829 //
ca58ed4e 1830
c9cbd2f2 1831 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
1832 if (!magF) AliError("Magneticd field - not initialized");
1833 Double_t bz = magF->SolenoidField(); //field in kGauss
1834 printf("bz: %lf\n",bz);
1835 AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
ca58ed4e 1836
c9cbd2f2 1837 TObjArray aTrk; // Original Track array of Aside
1838 TObjArray daTrk; // Distorted Track array of A side
1839 UShort_t *aId = new UShort_t[nTracks]; // A side Track ID
1840 TObjArray cTrk;
1841 TObjArray dcTrk;
1842 UShort_t *cId = new UShort_t [nTracks];
1843 Int_t id=0;
ca58ed4e 1844 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
7d85e147 1845 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
ca58ed4e 1846 fpt.SetParameters(7.24,0.120);
1847 fpt.SetNpx(10000);
1848 for(Int_t nt=0; nt<nTracks; nt++){
1849 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
7d85e147 1850 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
c9cbd2f2 1851 Double_t pt = fpt.GetRandom(); // momentum for f1
1852 // printf("phi %lf eta %lf pt %lf\n",phi,eta,pt);
ca58ed4e 1853 Short_t sign=1;
1854 if(gRandom->Rndm() < 0.5){
1855 sign =1;
1856 }else{
1857 sign=-1;
1858 }
1859
1860 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
1861 Double_t pxyz[3];
1862 pxyz[0]=pt*TMath::Cos(phi);
1863 pxyz[1]=pt*TMath::Sin(phi);
1864 pxyz[2]=pt*TMath::Tan(theta);
1865 Double_t cv[21]={0};
1866 AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
1867
1868 Double_t refX=1.;
1869 Int_t dir=-1;
1870 AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir, NULL);
1871 if (!td) continue;
1872 if (pcstream) (*pcstream)<<"track"<<
1873 "eta="<<eta<<
1874 "theta="<<theta<<
1875 "tOrig.="<<t<<
1876 "td.="<<td<<
1877 "\n";
7d85e147 1878 if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
ca58ed4e 1879 if (td){
c9cbd2f2 1880 daTrk.AddLast(td);
1881 aTrk.AddLast(t);
1882 Int_t nn=aTrk.GetEntriesFast();
1883 aId[nn]=id;
ca58ed4e 1884 }
7d85e147 1885 }else if(( eta<-0.07 )&&( eta>-etaCuts )){
ca58ed4e 1886 if (td){
c9cbd2f2 1887 dcTrk.AddLast(td);
1888 cTrk.AddLast(t);
1889 Int_t nn=cTrk.GetEntriesFast();
1890 cId[nn]=id;
ca58ed4e 1891 }
1892 }
c9cbd2f2 1893 id++;
ca58ed4e 1894 }// end of track loop
1895
1896 vertexer->SetTPCMode();
1897 vertexer->SetConstraintOff();
1898
c9cbd2f2 1899 aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));
1900 avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
1901 cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));
1902 cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
ca58ed4e 1903 if (pcstream) (*pcstream)<<"vertex"<<
1904 "x="<<orgVertex[0]<<
1905 "y="<<orgVertex[1]<<
1906 "z="<<orgVertex[2]<<
1907 "av.="<<&aV<< // distorted vertex A side
1908 "cv.="<<&cV<< // distroted vertex C side
1909 "avO.="<<&avOrg<< // original vertex A side
1910 "cvO.="<<&cvOrg<<
1911 "\n";
c9cbd2f2 1912 delete []aId;
1913 delete []cId;
ca58ed4e 1914}
f1817479 1915
1916void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
1917 //
1918 // make correction available for visualization using
1919 // TFormula, TFX and TTree::Draw
1920 // important in order to check corrections and also compute dervied variables
1921 // e.g correction partial derivatives
1922 //
1923 // NOTE - class is not owner of correction
1924 //
1925 if (!fgVisualCorrection) fgVisualCorrection=new TObjArray;
8847ede1 1926 if (position>=fgVisualCorrection->GetEntriesFast()) fgVisualCorrection->Expand(position*2);
f1817479 1927 fgVisualCorrection->AddAt(corr, position);
1928}
1929
1930
1931
287fbdfa 1932Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t r, Double_t kZ, Int_t axisType, Int_t corrType){
f1817479 1933 //
1934 // calculate the correction at given position - check the geffCorr
1935 //
1936 if (!fgVisualCorrection) return 0;
1937 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1938 if (!corr) return 0;
1939 Double_t phi=sector*TMath::Pi()/9.;
287fbdfa 1940 Double_t gx = r*TMath::Cos(phi);
1941 Double_t gy = r*TMath::Sin(phi);
1942 Double_t gz = r*kZ;
f1817479 1943 Int_t nsector=(gz>0) ? 0:18;
1944 //
1945 //
1946 //
1947 Float_t distPoint[3]={gx,gy,gz};
1948 corr->DistortPoint(distPoint, nsector);
1949 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1950 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1951 Double_t phi0=TMath::ATan2(gy,gx);
1952 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1953 if (axisType==0) return r1-r0;
1954 if (axisType==1) return (phi1-phi0)*r0;
1955 if (axisType==2) return distPoint[2]-gz;
1956 return phi1-phi0;
1957}
1958
1959Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
1960 //
1961 // return correction at given x,y,z
1962 //
1963 if (!fgVisualCorrection) return 0;
1964 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
1965 if (!corr) return 0;
1966 Double_t phi0= TMath::ATan2(gy,gx);
1967 Int_t nsector=(gz>0) ? 0:18;
1968 Float_t distPoint[3]={gx,gy,gz};
1969 corr->DistortPoint(distPoint, nsector);
1970 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
1971 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
1972 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
1973 if (axisType==0) return r1-r0;
1974 if (axisType==1) return (phi1-phi0)*r0;
1975 if (axisType==2) return distPoint[2]-gz;
1976 return phi1-phi0;
1977}