Restoring the generation of HTML documentation (Natalia) + related changes (PH)
[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
b4caed64 16// _________________________________________________________________
17//
18// Begin_Html
19// <h2> AliTPCCorrection class </h2>
20//
21// The AliTPCCorrection class provides a general framework to deal with space point distortions.
22// An correction class which inherits from here is for example AliTPCExBBShape or AliTPCExBTwist. <br>
23// General virtual functions are (for example) CorrectPoint(x,roc) where x is the vector of initial
24// positions in cartesian coordinates and roc represents the read-out chamber number according to
25// the offline numbering convention. The vector x is overwritten with the corrected coordinates. <br>
26// An alternative usage would be CorrectPoint(x,roc,dx), which leaves the vector x untouched, but
27// returns the distortions via the vector dx. <br>
28// This class is normally used via the general class AliTPCComposedCorrection.
29// <p>
30// Furthermore, the class contains basic geometrical descriptions like field cage radii
31// (fgkIFCRadius, fgkOFCRadius) and length (fgkTPCZ0) plus the voltages. Also, the definitions
32// of size and widths of the fulcrums building the grid of the final look-up table, which is
33// then interpolated, is defined in kNX and fgkXList).
34// <p>
35// All physics-model classes below are derived from this class in order to not duplicate code
36// and to allow a uniform treatment of all physics models.
37// <p>
38// <h3> Poisson solver </h3>
39// A numerical solver of the Poisson equation (relaxation technique) is implemented for 2-dimensional
40// geometries (r,z) as well as for 3-dimensional problems (r,$\phi$,z). The corresponding function
41// names are PoissonRelaxation?D. The relevant function arguments are the arrays of the boundary and
42// initial conditions (ArrayofArrayV, ArrayofChargeDensities) as well as the grid granularity which
43// is used during the calculation. These inputs can be chosen according to the needs of the physical
44// effect which is supposed to be simulated. In the 3D version, different symmetry conditions can be set
45// in order to reduce the calculation time (used in AliTPCFCVoltError3D).
46// <p>
47// <h3> Unified plotting functionality </h3>
48// Generic plot functions were implemented. They return a histogram pointer in the chosen plane of
49// the TPC drift volume with a selectable grid granularity and the magnitude of the correction vector.
50// For example, the function CreateHistoDZinXY(z,nx,ny) returns a 2-dimensional histogram which contains
51// the longitudinal corrections $dz$ in the (x,y)-plane at the given z position with the granularity of
52// nx and ny. The magnitude of the corrections is defined by the class from which this function is called.
53// In the same manner, standard plots for the (r,$\phi$)-plane and for the other corrections like $dr$ and $rd\phi$ are available
54// <p>
55// Note: This class is normally used via the class AliTPCComposedCorrection
56// End_Html
57//
6a1caa6b 58// Begin_Macro(source)
b4caed64 59// {
60// gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
6a1caa6b 61// TCanvas *c2 = new TCanvas("cAliTPCCorrection","cAliTPCCorrection",700,1050); c2->Divide(2,3);
b4caed64 62// AliTPCROCVoltError3D roc; // EXAMPLE PLOTS - SEE BELOW
63// roc.SetOmegaTauT1T2(0,1,1); // B=0
64// Float_t z0 = 1; // at +1 cm -> A side
65// c2->cd(1); roc.CreateHistoDRinXY(1.,300,300)->Draw("cont4z");
66// c2->cd(3);roc.CreateHistoDRPhiinXY(1.,300,300)->Draw("cont4z");
67// c2->cd(5);roc.CreateHistoDZinXY(1.,300,300)->Draw("cont4z");
68// Float_t phi0=0.5;
69// c2->cd(2);roc.CreateHistoDRinZR(phi0)->Draw("surf2");
70// c2->cd(4);roc.CreateHistoDRPhiinZR(phi0)->Draw("surf2");
71// c2->cd(6);roc.CreateHistoDZinZR(phi0)->Draw("surf2");
72// return c2;
73// }
74// End_Macro
75//
76// Begin_Html
77// <p>
78// Date: 27/04/2010 <br>
79// Authors: Magnus Mager, Stefan Rossegger, Jim Thomas
80// End_Html
81// _________________________________________________________________
82
83
be67055b 84#include "Riostream.h"
0116859c 85
86#include <TH2F.h>
87#include <TMath.h>
88#include <TROOT.h>
cf5b0aa0 89#include <TTreeStream.h>
ffab0c37 90#include <TTree.h>
91#include <TFile.h>
e527a1b9 92#include <TTimeStamp.h>
ffab0c37 93#include <AliCDBStorage.h>
94#include <AliCDBId.h>
95#include <AliCDBMetaData.h>
c9cbd2f2 96#include "TVectorD.h"
97#include "AliTPCParamSR.h"
7f4cb119 98
c9cbd2f2 99#include "AliTPCCorrection.h"
100#include "AliLog.h"
1b923461 101
1b923461 102#include "AliExternalTrackParam.h"
103#include "AliTrackPointArray.h"
104#include "TDatabasePDG.h"
105#include "AliTrackerBase.h"
106#include "AliTPCROC.h"
107#include "THnSparse.h"
108
c9cbd2f2 109#include "AliTPCLaserTrack.h"
110#include "AliESDVertex.h"
111#include "AliVertexerTracks.h"
112#include "TDatabasePDG.h"
113#include "TF1.h"
7f4cb119 114#include "TRandom.h"
c9cbd2f2 115
116#include "TDatabasePDG.h"
117
7f4cb119 118#include "AliTPCTransform.h"
119#include "AliTPCcalibDB.h"
120#include "AliTPCExB.h"
cf5b0aa0 121
c9cbd2f2 122#include "AliTPCRecoParam.h"
1b923461 123
0116859c 124
cf5b0aa0 125ClassImp(AliTPCCorrection)
126
f1817479 127
128TObjArray *AliTPCCorrection::fgVisualCorrection=0;
129// instance of correction for visualization
130
131
0116859c 132// FIXME: the following values should come from the database
c9cbd2f2 133const Double_t AliTPCCorrection::fgkTPCZ0 = 249.7; // nominal gating grid position
2b68ab9c 134const Double_t AliTPCCorrection::fgkIFCRadius= 83.5; // radius which renders the "18 rod manifold" best -> compare calc. of Jim Thomas
135// compare gkIFCRadius= 83.05: Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
c9cbd2f2 136const Double_t AliTPCCorrection::fgkOFCRadius= 254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
137const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
138const Double_t AliTPCCorrection::fgkCathodeV = -100000.0; // Cathode Voltage (volts)
139const Double_t AliTPCCorrection::fgkGG = -70.0; // Gating Grid voltage (volts)
0116859c 140
c9cbd2f2 141const Double_t AliTPCCorrection::fgkdvdE = 0.0024; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
0116859c 142
c9cbd2f2 143const Double_t AliTPCCorrection::fgkEM = -1.602176487e-19/9.10938215e-31; // charge/mass in [C/kg]
144const Double_t AliTPCCorrection::fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)]
c9cbd2f2 145
0116859c 146
147AliTPCCorrection::AliTPCCorrection()
c9cbd2f2 148 : TNamed("correction_unity","unity"),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 149{
150 //
151 // default constructor
152 //
f1817479 153 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
c9cbd2f2 154
35ae345f 155 InitLookUpfulcrums();
c9cbd2f2 156
0116859c 157}
158
159AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
c9cbd2f2 160: TNamed(name,title),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 161{
162 //
163 // default constructor, that set the name and title
164 //
f1817479 165 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
c9cbd2f2 166
35ae345f 167 InitLookUpfulcrums();
c9cbd2f2 168
0116859c 169}
170
171AliTPCCorrection::~AliTPCCorrection() {
172 //
173 // virtual destructor
174 //
175}
176
177void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
178 //
179 // Corrects the initial coordinates x (cartesian coordinates)
180 // according to the given effect (inherited classes)
181 // roc represents the TPC read out chamber (offline numbering convention)
182 //
183 Float_t dx[3];
184 GetCorrection(x,roc,dx);
185 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
186}
187
188void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
189 //
190 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
191 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
192 // roc represents the TPC read out chamber (offline numbering convention)
193 //
194 Float_t dx[3];
195 GetCorrection(x,roc,dx);
196 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
197}
198
199void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
200 //
201 // Distorts the initial coordinates x (cartesian coordinates)
202 // according to the given effect (inherited classes)
203 // roc represents the TPC read out chamber (offline numbering convention)
204 //
205 Float_t dx[3];
206 GetDistortion(x,roc,dx);
207 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
208}
209
46e89793 210void AliTPCCorrection::DistortPointLocal(Float_t x[],const Short_t roc) {
211 //
212 // Distorts the initial coordinates x (cartesian coordinates)
213 // according to the given effect (inherited classes)
214 // roc represents the TPC read out chamber (offline numbering convention)
215 //
216 Float_t gxyz[3]={0,0,0};
217 Double_t alpha = TMath::Pi()*(roc%18+0.5)/18;
218 Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
219 gxyz[0]= ca*x[0]+sa*x[1];
220 gxyz[1]= -sa*x[0]+ca*x[1];
221 gxyz[2]= x[2];
222 DistortPoint(gxyz,roc);
223 x[0]= ca*gxyz[0]-sa*gxyz[1];
224 x[1]= +sa*gxyz[0]+ca*gxyz[1];
225 x[2]= gxyz[2];
226}
227void AliTPCCorrection::CorrectPointLocal(Float_t x[],const Short_t roc) {
228 //
229 // Distorts the initial coordinates x (cartesian coordinates)
230 // according to the given effect (inherited classes)
231 // roc represents the TPC read out chamber (offline numbering convention)
232 //
233 Float_t gxyz[3]={0,0,0};
234 Double_t alpha = TMath::Pi()*(roc%18+0.5)/18;
235 Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
236 gxyz[0]= ca*x[0]+sa*x[1];
237 gxyz[1]= -sa*x[0]+ca*x[1];
238 gxyz[2]= x[2];
239 CorrectPoint(gxyz,roc);
240 x[0]= ca*gxyz[0]-sa*gxyz[1];
241 x[1]= sa*gxyz[0]+ca*gxyz[1];
242 x[2]= gxyz[2];
243}
244
0116859c 245void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
246 //
247 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
248 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
249 // roc represents the TPC read out chamber (offline numbering convention)
250 //
251 Float_t dx[3];
252 GetDistortion(x,roc,dx);
253 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
254}
255
256void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
257 //
258 // This function delivers the correction values dx in respect to the inital coordinates x
259 // roc represents the TPC read out chamber (offline numbering convention)
260 // Note: The dx is overwritten by the inherited effectice class ...
261 //
262 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
263}
264
265void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
266 //
267 // This function delivers the distortion values dx in respect to the inital coordinates x
268 // roc represents the TPC read out chamber (offline numbering convention)
269 //
270 GetCorrection(x,roc,dx);
271 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
272}
273
274void AliTPCCorrection::Init() {
275 //
276 // Initialization funtion (not used at the moment)
277 //
278}
279
e527a1b9 280void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
281 //
282 // Update function
283 //
284}
285
0116859c 286void AliTPCCorrection::Print(Option_t* /*option*/) const {
287 //
288 // Print function to check which correction classes are used
289 // option=="d" prints details regarding the setted magnitude
290 // option=="a" prints the C0 and C1 coefficents for calibration purposes
291 //
292 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
293}
294
534fd34a 295void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
0116859c 296 //
297 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
298 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
299 // calibration run
300 //
534fd34a 301 fT1=t1;
302 fT2=t2;
303 //SetOmegaTauT1T2(omegaTau, t1, t2);
0116859c 304}
305
306TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
307 //
308 // Simple plot functionality.
309 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
310 // in respect to position z within the XY plane.
311 // The histogramm has nx times ny entries.
312 //
c9cbd2f2 313 AliTPCParam* tpcparam = new AliTPCParamSR;
314
0116859c 315 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
316 nx,-250.,250.,ny,-250.,250.);
317 Float_t x[3],dx[3];
318 x[2]=z;
319 Int_t roc=z>0.?0:18; // FIXME
320 for (Int_t iy=1;iy<=ny;++iy) {
321 x[1]=h->GetYaxis()->GetBinCenter(iy);
322 for (Int_t ix=1;ix<=nx;++ix) {
323 x[0]=h->GetXaxis()->GetBinCenter(ix);
324 GetCorrection(x,roc,dx);
325 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
c9cbd2f2 326 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
0116859c 327 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
328 h->SetBinContent(ix,iy,r1-r0);
329 }
330 else
331 h->SetBinContent(ix,iy,0.);
332 }
333 }
c9cbd2f2 334 delete tpcparam;
0116859c 335 return h;
336}
337
338TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
339 //
340 // Simple plot functionality.
341 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
342 // in respect to position z within the XY plane.
343 // The histogramm has nx times ny entries.
344 //
345
c9cbd2f2 346 AliTPCParam* tpcparam = new AliTPCParamSR;
347
0116859c 348 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
349 nx,-250.,250.,ny,-250.,250.);
350 Float_t x[3],dx[3];
351 x[2]=z;
352 Int_t roc=z>0.?0:18; // FIXME
353 for (Int_t iy=1;iy<=ny;++iy) {
354 x[1]=h->GetYaxis()->GetBinCenter(iy);
355 for (Int_t ix=1;ix<=nx;++ix) {
356 x[0]=h->GetXaxis()->GetBinCenter(ix);
357 GetCorrection(x,roc,dx);
358 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
c9cbd2f2 359 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
0116859c 360 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
361 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
362
363 Float_t dphi=phi1-phi0;
364 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
365 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
366
367 h->SetBinContent(ix,iy,r0*dphi);
368 }
369 else
370 h->SetBinContent(ix,iy,0.);
371 }
372 }
c9cbd2f2 373 delete tpcparam;
374 return h;
375}
376
377TH2F* AliTPCCorrection::CreateHistoDZinXY(Float_t z,Int_t nx,Int_t ny) {
378 //
379 // Simple plot functionality.
380 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
381 // in respect to position z within the XY plane.
382 // The histogramm has nx times ny entries.
383 //
384
385 AliTPCParam* tpcparam = new AliTPCParamSR;
386
387 TH2F *h=CreateTH2F("dz_xy",GetTitle(),"x [cm]","y [cm]","dz [cm]",
388 nx,-250.,250.,ny,-250.,250.);
389 Float_t x[3],dx[3];
390 x[2]=z;
391 Int_t roc=z>0.?0:18; // FIXME
392 for (Int_t iy=1;iy<=ny;++iy) {
393 x[1]=h->GetYaxis()->GetBinCenter(iy);
394 for (Int_t ix=1;ix<=nx;++ix) {
395 x[0]=h->GetXaxis()->GetBinCenter(ix);
396 GetCorrection(x,roc,dx);
397 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
398 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
399 h->SetBinContent(ix,iy,dx[2]);
400 }
401 else
402 h->SetBinContent(ix,iy,0.);
403 }
404 }
405 delete tpcparam;
0116859c 406 return h;
407}
408
409TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
410 //
411 // Simple plot functionality.
412 // Returns a 2d hisogram which represents the corrections in r direction (dr)
413 // in respect to angle phi within the ZR plane.
414 // The histogramm has nx times ny entries.
415 //
416 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
417 nz,-250.,250.,nr,85.,250.);
418 Float_t x[3],dx[3];
419 for (Int_t ir=1;ir<=nr;++ir) {
420 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
421 x[0]=radius*TMath::Cos(phi);
422 x[1]=radius*TMath::Sin(phi);
423 for (Int_t iz=1;iz<=nz;++iz) {
424 x[2]=h->GetXaxis()->GetBinCenter(iz);
425 Int_t roc=x[2]>0.?0:18; // FIXME
426 GetCorrection(x,roc,dx);
427 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
428 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
429 h->SetBinContent(iz,ir,r1-r0);
430 }
431 }
0116859c 432 return h;
433
434}
435
436TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
437 //
438 // Simple plot functionality.
439 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
440 // in respect to angle phi within the ZR plane.
441 // The histogramm has nx times ny entries.
442 //
443 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
444 nz,-250.,250.,nr,85.,250.);
445 Float_t x[3],dx[3];
446 for (Int_t iz=1;iz<=nz;++iz) {
447 x[2]=h->GetXaxis()->GetBinCenter(iz);
448 Int_t roc=x[2]>0.?0:18; // FIXME
449 for (Int_t ir=1;ir<=nr;++ir) {
450 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
451 x[0]=radius*TMath::Cos(phi);
452 x[1]=radius*TMath::Sin(phi);
453 GetCorrection(x,roc,dx);
454 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
455 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
456 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
457
458 Float_t dphi=phi1-phi0;
459 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
460 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
461
462 h->SetBinContent(iz,ir,r0*dphi);
463 }
464 }
465 return h;
466}
467
c9cbd2f2 468TH2F* AliTPCCorrection::CreateHistoDZinZR(Float_t phi,Int_t nz,Int_t nr) {
469 //
470 // Simple plot functionality.
471 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
472 // in respect to angle phi within the ZR plane.
473 // The histogramm has nx times ny entries.
474 //
475 TH2F *h=CreateTH2F("dz_zr",GetTitle(),"z [cm]","r [cm]","dz [cm]",
476 nz,-250.,250.,nr,85.,250.);
477 Float_t x[3],dx[3];
478 for (Int_t ir=1;ir<=nr;++ir) {
479 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
480 x[0]=radius*TMath::Cos(phi);
481 x[1]=radius*TMath::Sin(phi);
482 for (Int_t iz=1;iz<=nz;++iz) {
483 x[2]=h->GetXaxis()->GetBinCenter(iz);
484 Int_t roc=x[2]>0.?0:18; // FIXME
485 GetCorrection(x,roc,dx);
486 h->SetBinContent(iz,ir,dx[2]);
487 }
488 }
489 return h;
490
491}
492
493
0116859c 494TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
495 const char *xlabel,const char *ylabel,const char *zlabel,
496 Int_t nbinsx,Double_t xlow,Double_t xup,
497 Int_t nbinsy,Double_t ylow,Double_t yup) {
498 //
499 // Helper function to create a 2d histogramm of given size
500 //
501
502 TString hname=name;
503 Int_t i=0;
504 if (gDirectory) {
505 while (gDirectory->FindObject(hname.Data())) {
506 hname =name;
507 hname+="_";
508 hname+=i;
509 ++i;
510 }
511 }
512 TH2F *h=new TH2F(hname.Data(),title,
513 nbinsx,xlow,xup,
514 nbinsy,ylow,yup);
515 h->GetXaxis()->SetTitle(xlabel);
516 h->GetYaxis()->SetTitle(ylabel);
517 h->GetZaxis()->SetTitle(zlabel);
518 h->SetStats(0);
519 return h;
520}
521
0116859c 522// Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
523
524void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
b1f0a2a5 525 const Double_t er[kNZ][kNR], Double_t &erValue ) {
0116859c 526 //
527 // Interpolate table - 2D interpolation
528 //
25732bff 529 Double_t saveEr[5] = {0,0,0,0,0};
0116859c 530
531 Search( kNZ, fgkZList, z, fJLow ) ;
532 Search( kNR, fgkRList, r, fKLow ) ;
533 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
534 if ( fKLow < 0 ) fKLow = 0 ;
535 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
536 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
537
538 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
b1f0a2a5 539 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
0116859c 540 }
b1f0a2a5 541 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
0116859c 542
543}
544
c9cbd2f2 545void AliTPCCorrection::Interpolate3DEdistortion( const Int_t order, const Double_t r, const Float_t phi, const Double_t z,
546 const Double_t er[kNZ][kNPhi][kNR], const Double_t ephi[kNZ][kNPhi][kNR], const Double_t ez[kNZ][kNPhi][kNR],
547 Double_t &erValue, Double_t &ephiValue, Double_t &ezValue) {
548 //
549 // Interpolate table - 3D interpolation
550 //
551
25732bff 552 Double_t saveEr[5]= {0,0,0,0,0};
553 Double_t savedEr[5]= {0,0,0,0,0} ;
554
555 Double_t saveEphi[5]= {0,0,0,0,0};
556 Double_t savedEphi[5]= {0,0,0,0,0} ;
557
558 Double_t saveEz[5]= {0,0,0,0,0};
559 Double_t savedEz[5]= {0,0,0,0,0} ;
c9cbd2f2 560
561 Search( kNZ, fgkZList, z, fILow ) ;
562 Search( kNPhi, fgkPhiList, z, fJLow ) ;
563 Search( kNR, fgkRList, r, fKLow ) ;
564
565 if ( fILow < 0 ) fILow = 0 ; // check if out of range
566 if ( fJLow < 0 ) fJLow = 0 ;
567 if ( fKLow < 0 ) fKLow = 0 ;
568
569 if ( fILow + order >= kNZ - 1 ) fILow = kNZ - 1 - order ;
570 if ( fJLow + order >= kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
571 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
572
573 for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
574 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
575 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r ) ;
576 saveEphi[j-fJLow] = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
577 saveEz[j-fJLow] = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r ) ;
578 }
579 savedEr[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi ) ;
580 savedEphi[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
581 savedEz[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi ) ;
582 }
583 erValue = Interpolate( &fgkZList[fILow], savedEr, order, z ) ;
584 ephiValue = Interpolate( &fgkZList[fILow], savedEphi, order, z ) ;
585 ezValue = Interpolate( &fgkZList[fILow], savedEz, order, z ) ;
586
587}
588
589Double_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
590 const Int_t nx, const Int_t ny, const Double_t xv[], const Double_t yv[],
591 const TMatrixD &array ) {
592 //
593 // Interpolate table (TMatrix format) - 2D interpolation
594 //
595
596 static Int_t jlow = 0, klow = 0 ;
25732bff 597 Double_t saveArray[5] = {0,0,0,0,0} ;
c9cbd2f2 598
599 Search( nx, xv, x, jlow ) ;
600 Search( ny, yv, y, klow ) ;
601 if ( jlow < 0 ) jlow = 0 ; // check if out of range
602 if ( klow < 0 ) klow = 0 ;
603 if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
604 if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
605
606 for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
607 {
608 Double_t *ajkl = &((TMatrixD&)array)(j,klow);
609 saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
610 }
611
612 return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
613
614}
615
616Double_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x, const Double_t y, const Double_t z,
617 const Int_t nx, const Int_t ny, const Int_t nz,
618 const Double_t xv[], const Double_t yv[], const Double_t zv[],
619 TMatrixD **arrayofArrays ) {
620 //
621 // Interpolate table (TMatrix format) - 3D interpolation
622 //
623
624 static Int_t ilow = 0, jlow = 0, klow = 0 ;
25732bff 625 Double_t saveArray[5]= {0,0,0,0,0};
626 Double_t savedArray[5]= {0,0,0,0,0} ;
c9cbd2f2 627
628 Search( nx, xv, x, ilow ) ;
629 Search( ny, yv, y, jlow ) ;
630 Search( nz, zv, z, klow ) ;
631
632 if ( ilow < 0 ) ilow = 0 ; // check if out of range
633 if ( jlow < 0 ) jlow = 0 ;
634 if ( klow < 0 ) klow = 0 ;
635
636 if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
637 if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
638 if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
639
640 for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
641 {
642 TMatrixD &table = *arrayofArrays[k] ;
643 for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
644 {
645 saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
646 }
647 savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
648 }
649 return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
650
651}
652
0116859c 653Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
b1f0a2a5 654 const Int_t order, const Double_t x ) {
0116859c 655 //
656 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
657 //
658
659 Double_t y ;
660 if ( order == 2 ) { // Quadratic Interpolation = 2
661 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
662 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
663 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
664 } else { // Linear Interpolation = 1
665 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
666 }
667
668 return (y);
669
670}
671
2bf29b72 672Float_t AliTPCCorrection::Interpolate2DTable( const Int_t order, const Double_t x, const Double_t y,
673 const Int_t nx, const Int_t ny, const Double_t xv[], const Double_t yv[],
674 const TMatrixF &array ) {
675 //
676 // Interpolate table (TMatrix format) - 2D interpolation
677 // Float version (in order to decrease the OCDB size)
678 //
679
680 static Int_t jlow = 0, klow = 0 ;
681 Float_t saveArray[5] = {0.,0.,0.,0.,0.} ;
682
683 Search( nx, xv, x, jlow ) ;
684 Search( ny, yv, y, klow ) ;
685 if ( jlow < 0 ) jlow = 0 ; // check if out of range
686 if ( klow < 0 ) klow = 0 ;
687 if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
688 if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
689
690 for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
691 {
692 Float_t *ajkl = &((TMatrixF&)array)(j,klow);
693 saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
694 }
695
696 return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
697
698}
699
700Float_t AliTPCCorrection::Interpolate3DTable( const Int_t order, const Double_t x, const Double_t y, const Double_t z,
701 const Int_t nx, const Int_t ny, const Int_t nz,
702 const Double_t xv[], const Double_t yv[], const Double_t zv[],
703 TMatrixF **arrayofArrays ) {
704 //
705 // Interpolate table (TMatrix format) - 3D interpolation
706 // Float version (in order to decrease the OCDB size)
707 //
708
709 static Int_t ilow = 0, jlow = 0, klow = 0 ;
710 Float_t saveArray[5]= {0.,0.,0.,0.,0.};
711 Float_t savedArray[5]= {0.,0.,0.,0.,0.} ;
712
713 Search( nx, xv, x, ilow ) ;
714 Search( ny, yv, y, jlow ) ;
715 Search( nz, zv, z, klow ) ;
716
717 if ( ilow < 0 ) ilow = 0 ; // check if out of range
718 if ( jlow < 0 ) jlow = 0 ;
719 if ( klow < 0 ) klow = 0 ;
720
721 if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
722 if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
723 if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
724
725 for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
726 {
727 TMatrixF &table = *arrayofArrays[k] ;
728 for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
729 {
730 saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
731 }
732 savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
733 }
734 return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
735
736}
737Float_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Float_t yArray[],
738 const Int_t order, const Double_t x ) {
739 //
740 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
741 // Float version (in order to decrease the OCDB size)
742 //
743
744 Float_t y ;
745 if ( order == 2 ) { // Quadratic Interpolation = 2
746 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
747 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
748 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
749 } else { // Linear Interpolation = 1
750 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
751 }
752
753 return (y);
754
755}
756
757
0116859c 758
b1f0a2a5 759void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
0116859c 760 //
761 // Search an ordered table by starting at the most recently used point
762 //
763
764 Long_t middle, high ;
765 Int_t ascend = 0, increment = 1 ;
766
767 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
768
769 if ( low < 0 || low > n-1 ) {
770 low = -1 ; high = n ;
771 } else { // Ordered Search phase
772 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
773 if ( low == n-1 ) return ;
774 high = low + 1 ;
775 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
776 low = high ;
777 increment *= 2 ;
778 high = low + increment ;
779 if ( high > n-1 ) { high = n ; break ; }
780 }
781 } else {
782 if ( low == 0 ) { low = -1 ; return ; }
783 high = low - 1 ;
784 while ( (Int_t)( x < xArray[low] ) == ascend ) {
785 high = low ;
786 increment *= 2 ;
787 if ( increment >= high ) { low = -1 ; break ; }
788 else low = high - increment ;
789 }
790 }
791 }
792
793 while ( (high-low) != 1 ) { // Binary Search Phase
794 middle = ( high + low ) / 2 ;
795 if ( (Int_t)( x >= xArray[middle] ) == ascend )
796 low = middle ;
797 else
798 high = middle ;
799 }
800
801 if ( x == xArray[n-1] ) low = n-2 ;
802 if ( x == xArray[0] ) low = 0 ;
803
804}
805
35ae345f 806void AliTPCCorrection::InitLookUpfulcrums() {
807 //
808 // Initialization of interpolation points - for main look up table
809 // (course grid in the middle, fine grid on the borders)
810 //
811
812 AliTPCROC * roc = AliTPCROC::Instance();
813 const Double_t rLow = TMath::Floor(roc->GetPadRowRadii(0,0))-1; // first padRow plus some margin
814
815 // fulcrums in R
816 fgkRList[0] = rLow;
817 for (Int_t i = 1; i<kNR; i++) {
818 fgkRList[i] = fgkRList[i-1] + 3.5; // 3.5 cm spacing
819 if (fgkRList[i]<90 ||fgkRList[i]>245)
820 fgkRList[i] = fgkRList[i-1] + 0.5; // 0.5 cm spacing
821 else if (fgkRList[i]<100 || fgkRList[i]>235)
822 fgkRList[i] = fgkRList[i-1] + 1.5; // 1.5 cm spacing
823 else if (fgkRList[i]<120 || fgkRList[i]>225)
824 fgkRList[i] = fgkRList[i-1] + 2.5; // 2.5 cm spacing
825 }
826
827 // fulcrums in Z
828 fgkZList[0] = -249.5;
829 fgkZList[kNZ-1] = 249.5;
830 for (Int_t j = 1; j<kNZ/2; j++) {
831 fgkZList[j] = fgkZList[j-1];
832 if (TMath::Abs(fgkZList[j])< 0.15)
833 fgkZList[j] = fgkZList[j-1] + 0.09; // 0.09 cm spacing
834 else if(TMath::Abs(fgkZList[j])< 0.6)
835 fgkZList[j] = fgkZList[j-1] + 0.4; // 0.4 cm spacing
836 else if (TMath::Abs(fgkZList[j])< 2.5 || TMath::Abs(fgkZList[j])>248)
837 fgkZList[j] = fgkZList[j-1] + 0.5; // 0.5 cm spacing
838 else if (TMath::Abs(fgkZList[j])<10 || TMath::Abs(fgkZList[j])>235)
839 fgkZList[j] = fgkZList[j-1] + 1.5; // 1.5 cm spacing
840 else if (TMath::Abs(fgkZList[j])<25 || TMath::Abs(fgkZList[j])>225)
841 fgkZList[j] = fgkZList[j-1] + 2.5; // 2.5 cm spacing
842 else
843 fgkZList[j] = fgkZList[j-1] + 4; // 4 cm spacing
844
845 fgkZList[kNZ-j-1] = -fgkZList[j];
846 }
847
848 // fulcrums in phi
849 for (Int_t k = 0; k<kNPhi; k++)
850 fgkPhiList[k] = TMath::TwoPi()*k/(kNPhi-1);
851
852
853}
854
855
c9cbd2f2 856void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
857 TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
858 const Int_t rows, const Int_t columns, const Int_t iterations,
859 const Bool_t rocDisplacement ) {
1b923461 860 //
861 // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
862 //
863 // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
864 // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
865 // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
866 // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
867 // you wish to solve Laplaces equation however it should not contain random numbers or you will get
868 // random numbers back as a solution.
869 // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
870 // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
871 // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
872 // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
873 // number of points and solves the problem again. This happens several times until the maximum number
874 // of points has been included in the array.
875 //
876 // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
877 // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
878 //
c9cbd2f2 879 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
880 //
1b923461 881 // Original code by Jim Thomas (STAR TPC Collaboration)
882 //
883
884 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
885
886 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
887 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
888 const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
889
890 TMatrixD arrayEr(rows,columns) ;
891 TMatrixD arrayEz(rows,columns) ;
892
893 //Check that number of rows and columns is suitable for a binary expansion
894
895 if ( !IsPowerOfTwo(rows-1) ) {
896 AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
897 return;
898 }
899 if ( !IsPowerOfTwo(columns-1) ) {
900 AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
901 return;
902 }
903
904 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
905 // Allow for different size grid spacing in R and Z directions
906 // Use a binary expansion of the size of the matrix to speed up the solution of the problem
907
908 Int_t iOne = (rows-1)/4 ;
909 Int_t jOne = (columns-1)/4 ;
910 // Solve for N in 2**N, add one.
911 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
912
913 for ( Int_t count = 0 ; count < loops ; count++ ) {
914 // Loop while the matrix expands & the resolution increases.
915
916 Float_t tempGridSizeR = gridSizeR * iOne ;
917 Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
918 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
919
920 // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
921 std::vector<float> coef1(rows) ;
922 std::vector<float> coef2(rows) ;
923
924 for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
925 Float_t radius = fgkIFCRadius + i*gridSizeR ;
926 coef1[i] = 1.0 + tempGridSizeR/(2*radius);
927 coef2[i] = 1.0 - tempGridSizeR/(2*radius);
928 }
929
930 TMatrixD sumChargeDensity(rows,columns) ;
931
932 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
933 Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
934 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
935 if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
936 else {
937 // Add up all enclosed charge density contributions within 1/2 unit in all directions
938 Float_t weight = 0.0 ;
939 Float_t sum = 0.0 ;
940 sumChargeDensity(i,j) = 0.0 ;
941 for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
942 for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
943 if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
944 else
945 weight = 1.0 ;
946 // Note that this is cylindrical geometry
947 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
948 sum += weight*radius ;
949 }
950 }
951 sumChargeDensity(i,j) /= sum ;
952 }
953 sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
954 }
955 }
956
957 for ( Int_t k = 1 ; k <= iterations; k++ ) {
958 // Solve Poisson's Equation
959 // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
960 // as interations increase.
961 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
962 Float_t overRelaxM1 = overRelax - 1.0 ;
963 Float_t overRelaxtempFourth, overRelaxcoef5 ;
964 overRelaxtempFourth = overRelax * tempFourth ;
965 overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
966
967 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
968 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
969
970 arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
971 + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
972 - overRelaxcoef5 * arrayV(i,j)
973 + coef1[i] * arrayV(i+iOne,j)
974 + sumChargeDensity(i,j)
975 ) * overRelaxtempFourth;
976 }
977 }
978
979 if ( k == iterations ) {
980 // After full solution is achieved, copy low resolution solution into higher res array
981 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
982 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
983
984 if ( iOne > 1 ) {
985 arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
986 if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
987 }
988 if ( jOne > 1 ) {
989 arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
990 if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
991 }
992 if ( iOne > 1 && jOne > 1 ) {
993 arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
994 if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
995 if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
996 // Note that this leaves a point at the upper left and lower right corners uninitialized.
997 // -> Not a big deal.
998 }
999
1000 }
1001 }
1002 }
1003
1004 }
1005
1006 iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
1007 jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
1008
c9cbd2f2 1009 sumChargeDensity.Clear();
1b923461 1010 }
1011
1012 // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
1013 for ( Int_t j = 0 ; j < columns ; j++ ) {
1014 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1015 arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1016 arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1017 }
1018
1019 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1020 for ( Int_t i = 0 ; i < rows ; i++) {
1021 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1022 arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1023 arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1024 }
1025
1026 for ( Int_t i = 0 ; i < rows ; i++) {
1027 // Note: go back and compare to old version of this code. See notes below.
1028 // JT Test ... attempt to divide by real Ez not Ez to first order
1029 for ( Int_t j = 0 ; j < columns ; j++ ) {
1030 arrayEz(i,j) += ezField;
1031 // This adds back the overall Z gradient of the field (main E field component)
1032 }
1033 // Warning: (-=) assumes you are using an error potetial without the overall Field included
1034 }
1035
1036 // Integrate Er/Ez from Z to zero
1037 for ( Int_t j = 0 ; j < columns ; j++ ) {
1038 for ( Int_t i = 0 ; i < rows ; i++ ) {
c9cbd2f2 1039
1b923461 1040 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1041 arrayErOverEz(i,j) = 0.0 ;
c9cbd2f2 1042 arrayDeltaEz(i,j) = 0.0 ;
1043
1b923461 1044 for ( Int_t k = j ; k < columns ; k++ ) {
1045 arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
c9cbd2f2 1046 arrayDeltaEz(i,j) += index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
1b923461 1047 if ( index != 4 ) index = 4; else index = 2 ;
1048 }
c9cbd2f2 1049 if ( index == 4 ) {
1050 arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
1051 arrayDeltaEz(i,j) -= (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
1052 }
1053 if ( index == 2 ) {
1054 arrayErOverEz(i,j) += (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
1055 -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
1056 arrayDeltaEz(i,j) += (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
1057 -2.5*(arrayEz(i,columns-1)-ezField));
1058 }
1059 if ( j == columns-2 ) {
1060 arrayErOverEz(i,j) = (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
1061 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
1062 arrayDeltaEz(i,j) = (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
1063 +1.5*(arrayEz(i,columns-1)-ezField) ) ;
1064 }
1065 if ( j == columns-1 ) {
1066 arrayErOverEz(i,j) = 0.0 ;
1067 arrayDeltaEz(i,j) = 0.0 ;
1068 }
1b923461 1069 }
1070 }
1071
c9cbd2f2 1072 // calculate z distortion from the integrated Delta Ez residuals
1073 // and include the aquivalence (Volt to cm) of the ROC shift !!
1074
1075 for ( Int_t j = 0 ; j < columns ; j++ ) {
1076 for ( Int_t i = 0 ; i < rows ; i++ ) {
1077
1078 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1079 arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
1080
1081 // ROC Potential in cm aquivalent
1082 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
1083 if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift; // add the ROC misaligment
1084
1085 }
1086 }
1087
1088 arrayEr.Clear();
1089 arrayEz.Clear();
1090
1b923461 1091}
1092
c9cbd2f2 1093void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
1094 TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
1095 const Int_t rows, const Int_t columns, const Int_t phislices,
1096 const Float_t deltaphi, const Int_t iterations, const Int_t symmetry,
1097 Bool_t rocDisplacement ) {
1098 //
1099 // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
1100 //
1101 // NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
1102 // The number of rows and COLUMNS can be different.
1103 //
1104 // ROWS == 2**M + 1
1105 // COLUMNS == 2**N + 1
1106 // PHISLICES == Arbitrary but greater than 3
1107 //
1108 // DeltaPhi in Radians
1109 //
1110 // SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
1111 // = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
1112 //
1113 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
1114
1115 const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
1116
1117 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
1118 const Float_t gridSizePhi = deltaphi ;
1119 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
1120 const Float_t ratioPhi = gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
1121 const Float_t ratioZ = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
1122
1123 TMatrixD arrayE(rows,columns) ;
1124
1125 // Check that the number of rows and columns is suitable for a binary expansion
1126 if ( !IsPowerOfTwo((rows-1)) ) {
1127 AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
1128 return; }
1129 if ( !IsPowerOfTwo((columns-1)) ) {
1130 AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
1131 return; }
1132 if ( phislices <= 3 ) {
1133 AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
1134 return; }
1135 if ( phislices > 1000 ) {
1136 AliError("Poisson3D phislices > 1000 is not allowed (nor wise) ");
1137 return; }
1138
1139 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
1140 // Allow for different size grid spacing in R and Z directions
1141 // Use a binary expansion of the matrix to speed up the solution of the problem
1142
1143 Int_t loops, mplus, mminus, signplus, signminus ;
1144 Int_t ione = (rows-1)/4 ;
1145 Int_t jone = (columns-1)/4 ;
1146 loops = TMath::Max(ione, jone) ; // Calculate the number of loops for the binary expansion
1147 loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ; // Solve for N in 2**N
1148
1149 TMatrixD* arrayofSumChargeDensities[1000] ; // Create temporary arrays to store low resolution charge arrays
1150
1151 for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
1152
1153 for ( Int_t count = 0 ; count < loops ; count++ ) { // START the master loop and do the binary expansion
1154
1155 Float_t tempgridSizeR = gridSizeR * ione ;
1156 Float_t tempratioPhi = ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
1157 Float_t tempratioZ = ratioZ * ione * ione / ( jone * jone ) ;
1158
1159 std::vector<float> coef1(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1160 std::vector<float> coef2(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1161 std::vector<float> coef3(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1162 std::vector<float> coef4(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1163
1164 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1165 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1166 coef1[i] = 1.0 + tempgridSizeR/(2*radius);
1167 coef2[i] = 1.0 - tempgridSizeR/(2*radius);
1168 coef3[i] = tempratioPhi/(radius*radius);
1169 coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
1170 }
1171
1172 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1173 TMatrixD &chargeDensity = *arrayofChargeDensities[m] ;
1174 TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
1175 for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
1176 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1177 for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
1178 if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1179 else { // Add up all enclosed charge density contributions within 1/2 unit in all directions
1180 Float_t weight = 0.0 ;
1181 Float_t sum = 0.0 ;
1182 sumChargeDensity(i,j) = 0.0 ;
1183 for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
1184 for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
1185 if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
1186 else
1187 weight = 1.0 ;
1188 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1189 sum += weight*radius ;
1190 }
1191 }
1192 sumChargeDensity(i,j) /= sum ;
1193 }
1194 sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
1195 }
1196 }
1197 }
1198
1199 for ( Int_t k = 1 ; k <= iterations; k++ ) {
1200
1201 // over-relaxation index, >= 1 but < 2
1202 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1203 Float_t overRelaxM1 = overRelax - 1.0 ;
1204
1205 std::vector<float> overRelaxcoef4(rows) ; // Do this the standard C++ way to avoid gcc extensions
1206 std::vector<float> overRelaxcoef5(rows) ; // Do this the standard C++ way to avoid gcc extensions
1207
1208 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1209 overRelaxcoef4[i] = overRelax * coef4[i] ;
1210 overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
1211 }
1212
1213 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1214
1215 mplus = m + 1; signplus = 1 ;
1216 mminus = m - 1 ; signminus = 1 ;
1217 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1218 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1219 if ( mminus < 0 ) mminus = 1 ;
1220 }
1221 else if (symmetry==-1) { // Anti-symmetry in phi
1222 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1223 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1224 }
1225 else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1226 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1227 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1228 }
1229 TMatrixD& arrayV = *arrayofArrayV[m] ;
1230 TMatrixD& arrayVP = *arrayofArrayV[mplus] ;
1231 TMatrixD& arrayVM = *arrayofArrayV[mminus] ;
1232 TMatrixD& sumChargeDensity = *arrayofSumChargeDensities[m] ;
1233
1234 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1235 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1236
1237 arrayV(i,j) = ( coef2[i] * arrayV(i-ione,j)
1238 + tempratioZ * ( arrayV(i,j-jone) + arrayV(i,j+jone) )
1239 - overRelaxcoef5[i] * arrayV(i,j)
1240 + coef1[i] * arrayV(i+ione,j)
1241 + coef3[i] * ( signplus*arrayVP(i,j) + signminus*arrayVM(i,j) )
1242 + sumChargeDensity(i,j)
1243 ) * overRelaxcoef4[i] ;
1244 // Note: over-relax the solution at each step. This speeds up the convergance.
1245
1246 }
1247 }
1248
1249 if ( k == iterations ) { // After full solution is achieved, copy low resolution solution into higher res array
1250 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1251 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1252
1253 if ( ione > 1 ) {
1254 arrayV(i+ione/2,j) = ( arrayV(i+ione,j) + arrayV(i,j) ) / 2 ;
1255 if ( i == ione ) arrayV(i-ione/2,j) = ( arrayV(0,j) + arrayV(ione,j) ) / 2 ;
1256 }
1257 if ( jone > 1 ) {
1258 arrayV(i,j+jone/2) = ( arrayV(i,j+jone) + arrayV(i,j) ) / 2 ;
1259 if ( j == jone ) arrayV(i,j-jone/2) = ( arrayV(i,0) + arrayV(i,jone) ) / 2 ;
1260 }
1261 if ( ione > 1 && jone > 1 ) {
1262 arrayV(i+ione/2,j+jone/2) = ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
1263 if ( i == ione ) arrayV(i-ione/2,j-jone/2) = ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
1264 if ( j == jone ) arrayV(i-ione/2,j-jone/2) = ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
1265 // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
1266 }
1267 }
1268 }
1269 }
1270
1271 }
1272 }
1273
1274 ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
1275 jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
1276
1277 }
1278
1279 //Differentiate V(r) and solve for E(r) using special equations for the first and last row
1280 //Integrate E(r)/E(z) from point of origin to pad plane
1281
1282 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1283 TMatrixD& arrayV = *arrayofArrayV[m] ;
1284 TMatrixD& eroverEz = *arrayofEroverEz[m] ;
1285
1286 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1287
1288 // Differentiate in R
1289 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1290 arrayE(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1291 arrayE(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1292 // Integrate over Z
1293 for ( Int_t i = 0 ; i < rows ; i++ ) {
1294 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1295 eroverEz(i,j) = 0.0 ;
1296 for ( Int_t k = j ; k < columns ; k++ ) {
1297
1298 eroverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1299 if ( index != 4 ) index = 4; else index = 2 ;
1300 }
1301 if ( index == 4 ) eroverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1302 if ( index == 2 ) eroverEz(i,j) +=
1303 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1304 if ( j == columns-2 ) eroverEz(i,j) =
1305 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1306 if ( j == columns-1 ) eroverEz(i,j) = 0.0 ;
1307 }
1308 }
1309 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1310 // eroverEz.Draw("surf") ; } // JT test
1311 }
1312
1313 //Differentiate V(r) and solve for E(phi)
1314 //Integrate E(phi)/E(z) from point of origin to pad plane
1315
1316 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1317
1318 mplus = m + 1; signplus = 1 ;
1319 mminus = m - 1 ; signminus = 1 ;
1320 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1321 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1322 if ( mminus < 0 ) mminus = 1 ;
1323 }
1324 else if (symmetry==-1) { // Anti-symmetry in phi
1325 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1326 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1327 }
1328 else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
1329 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1330 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1331 }
1332 TMatrixD &arrayVP = *arrayofArrayV[mplus] ;
1333 TMatrixD &arrayVM = *arrayofArrayV[mminus] ;
1334 TMatrixD &ePhioverEz = *arrayofEPhioverEz[m] ;
1335 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1336 // Differentiate in Phi
1337 for ( Int_t i = 0 ; i < rows ; i++ ) {
1338 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1339 arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
1340 }
1341 // Integrate over Z
1342 for ( Int_t i = 0 ; i < rows ; i++ ) {
1343 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1344 ePhioverEz(i,j) = 0.0 ;
1345 for ( Int_t k = j ; k < columns ; k++ ) {
1346
1347 ePhioverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1348 if ( index != 4 ) index = 4; else index = 2 ;
1349 }
1350 if ( index == 4 ) ePhioverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1351 if ( index == 2 ) ePhioverEz(i,j) +=
1352 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1353 if ( j == columns-2 ) ePhioverEz(i,j) =
1354 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1355 if ( j == columns-1 ) ePhioverEz(i,j) = 0.0 ;
1356 }
1357 }
1358 // if ( m == 5 ) { TCanvas* c2 = new TCanvas("arrayE","arrayE",50,50,840,600) ; c2 -> cd() ;
1359 // arrayE.Draw("surf") ; } // JT test
1360 }
1361
1362
1363 // Differentiate V(r) and solve for E(z) using special equations for the first and last row
1364 // Integrate (E(z)-Ezstd) from point of origin to pad plane
1365
1366 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1367 TMatrixD& arrayV = *arrayofArrayV[m] ;
1368 TMatrixD& deltaEz = *arrayofDeltaEz[m] ;
1369
1370 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1371 for ( Int_t i = 0 ; i < rows ; i++) {
1372 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1373 arrayE(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1374 arrayE(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1375 }
1376
1377 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1378 // Integrate over Z
1379 for ( Int_t i = 0 ; i < rows ; i++ ) {
1380 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1381 deltaEz(i,j) = 0.0 ;
1382 for ( Int_t k = j ; k < columns ; k++ ) {
1383 deltaEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k) ;
1384 if ( index != 4 ) index = 4; else index = 2 ;
1385 }
1386 if ( index == 4 ) deltaEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1) ;
1387 if ( index == 2 ) deltaEz(i,j) +=
1388 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
1389 if ( j == columns-2 ) deltaEz(i,j) =
1390 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
1391 if ( j == columns-1 ) deltaEz(i,j) = 0.0 ;
1392 }
1393 }
1394 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1395 // eroverEz.Draw("surf") ; } // JT test
1396
1397 // calculate z distortion from the integrated Delta Ez residuals
1398 // and include the aquivalence (Volt to cm) of the ROC shift !!
1399
1400 for ( Int_t j = 0 ; j < columns ; j++ ) {
1401 for ( Int_t i = 0 ; i < rows ; i++ ) {
1402
1403 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1404 deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
1405
1406 // ROC Potential in cm aquivalent
1407 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
1408 if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift; // add the ROC misaligment
1409
1410 }
1411 }
1412
1413 } // end loop over phi
1414
1415
1416
1417 for ( Int_t k = 0 ; k < phislices ; k++ )
1418 {
1419 arrayofSumChargeDensities[k]->Delete() ;
1420 }
1421
1422
1423
1424 arrayE.Clear();
1425}
1b923461 1426
1427
710bda39 1428Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1b923461 1429 //
1430 // Helperfunction: Check if integer is a power of 2
1431 //
1432 Int_t j = 0;
1433 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
1434 if ( j == 1 ) return(1) ; // True
1435 return(0) ; // False
1436}
1437
cf5b0aa0 1438
b1f0a2a5 1439AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
cf5b0aa0 1440 //
1441 // Fit the track parameters - without and with distortion
1442 // 1. Space points in the TPC are simulated along the trajectory
1443 // 2. Space points distorted
1444 // 3. Fits the non distorted and distroted track to the reference plane at refX
1445 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
1446 //
1447 // trackIn - input track parameters
1448 // refX - reference X to fit the track
1449 // dir - direction - out=1 or in=-1
1450 // pcstream - debug streamer to check the results
1451 //
cad404e1 1452 // see AliExternalTrackParam.h documentation:
1453 // track1.fP[0] - local y (rphi)
1454 // track1.fP[1] - z
1455 // track1.fP[2] - sinus of local inclination angle
1456 // track1.fP[3] - tangent of deep angle
1457 // track1.fP[4] - 1/pt
1b923461 1458
cf5b0aa0 1459 AliTPCROC * roc = AliTPCROC::Instance();
1460 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
1461 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
1462 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
cf5b0aa0 1463 const Double_t kMaxSnp = 0.85;
1464 const Double_t kSigmaY=0.1;
1465 const Double_t kSigmaZ=0.1;
ca58ed4e 1466 const Double_t kMaxR=500;
1467 const Double_t kMaxZ=500;
46e89793 1468
cfe2c39a 1469 const Double_t kMaxZ0=220;
1470 const Double_t kZcut=3;
cf5b0aa0 1471 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
ca58ed4e 1472 Int_t npoints1=0;
1473 Int_t npoints2=0;
cf5b0aa0 1474
be67055b 1475 AliExternalTrackParam track(trackIn); //
cf5b0aa0 1476 // generate points
1477 AliTrackPointArray pointArray0(npoints0);
1478 AliTrackPointArray pointArray1(npoints0);
1479 Double_t xyz[3];
cfe2c39a 1480 if (!AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,5,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 1481 //
1482 // simulate the track
1483 Int_t npoints=0;
1484 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
1485 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
cfe2c39a 1486 if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) return 0;
cf5b0aa0 1487 track.GetXYZ(xyz);
cfe2c39a 1488 xyz[0]+=gRandom->Gaus(0,0.000005);
1489 xyz[1]+=gRandom->Gaus(0,0.000005);
1490 xyz[2]+=gRandom->Gaus(0,0.000005);
1491 if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
46e89793 1492 if (TMath::Abs(track.GetX())<kRTPC0) continue;
1493 if (TMath::Abs(track.GetX())>kRTPC1) continue;
cf5b0aa0 1494 AliTrackPoint pIn0; // space point
1495 AliTrackPoint pIn1;
ffab0c37 1496 Int_t sector= (xyz[2]>0)? 0:18;
cf5b0aa0 1497 pointArray0.GetPoint(pIn0,npoints);
1498 pointArray1.GetPoint(pIn1,npoints);
1499 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
1500 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
ffab0c37 1501 DistortPoint(distPoint, sector);
cf5b0aa0 1502 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
1503 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
1504 //
1505 track.Rotate(alpha);
1506 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
1507 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
1508 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
1509 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
1510 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
1511 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
1512 pointArray0.AddPoint(npoints, &pIn0);
1513 pointArray1.AddPoint(npoints, &pIn1);
1514 npoints++;
1515 if (npoints>=npoints0) break;
1516 }
cfe2c39a 1517 if (npoints<npoints0/4.) return 0;
cf5b0aa0 1518 //
1519 // refit track
1520 //
1521 AliExternalTrackParam *track0=0;
1522 AliExternalTrackParam *track1=0;
1523 AliTrackPoint point1,point2,point3;
1524 if (dir==1) { //make seed inner
1525 pointArray0.GetPoint(point1,1);
cfe2c39a 1526 pointArray0.GetPoint(point2,11);
1527 pointArray0.GetPoint(point3,21);
cf5b0aa0 1528 }
1529 if (dir==-1){ //make seed outer
cfe2c39a 1530 pointArray0.GetPoint(point1,npoints-21);
1531 pointArray0.GetPoint(point2,npoints-11);
cf5b0aa0 1532 pointArray0.GetPoint(point3,npoints-1);
46e89793 1533 }
1534 if ((TMath::Abs(point1.GetX()-point3.GetX())+TMath::Abs(point1.GetY()-point3.GetY()))<10){
1535 printf("fit points not properly initialized\n");
1536 return 0;
1537 }
cf5b0aa0 1538 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
1539 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
cfe2c39a 1540 track0->ResetCovariance(10);
1541 track1->ResetCovariance(10);
1542 if (TMath::Abs(AliTrackerBase::GetBz())<0.01){
1543 ((Double_t*)track0->GetParameter())[4]= trackIn.GetParameter()[4];
1544 ((Double_t*)track1->GetParameter())[4]= trackIn.GetParameter()[4];
1545 }
cf5b0aa0 1546 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
8b63d99c 1547 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
cf5b0aa0 1548 //
1549 AliTrackPoint pIn0;
1550 AliTrackPoint pIn1;
1551 pointArray0.GetPoint(pIn0,ipoint);
1552 pointArray1.GetPoint(pIn1,ipoint);
1553 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
1554 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
46e89793 1555 if (TMath::Abs(prot0.GetX())<kRTPC0) continue;
1556 if (TMath::Abs(prot0.GetX())>kRTPC1) continue;
cf5b0aa0 1557 //
cfe2c39a 1558 if (!AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
1559 if (!AliTrackerBase::PropagateTrackTo(track1,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
ca58ed4e 1560 if (TMath::Abs(track0->GetZ())>kMaxZ) break;
1561 if (TMath::Abs(track0->GetX())>kMaxR) break;
1562 if (TMath::Abs(track1->GetZ())>kMaxZ) break;
1563 if (TMath::Abs(track1->GetX())>kMaxR) break;
cfe2c39a 1564 if (dir>0 && track1->GetX()>refX) continue;
1565 if (dir<0 && track1->GetX()<refX) continue;
1566 if (TMath::Abs(track1->GetZ())<kZcut)continue;
8b63d99c 1567 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
cf5b0aa0 1568 //
1569 Double_t pointPos[2]={0,0};
1570 Double_t pointCov[3]={0,0,0};
1571 pointPos[0]=prot0.GetY();//local y
1572 pointPos[1]=prot0.GetZ();//local z
1573 pointCov[0]=prot0.GetCov()[3];//simay^2
1574 pointCov[1]=prot0.GetCov()[4];//sigmayz
1575 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
ca58ed4e 1576 if (!track0->Update(pointPos,pointCov)) break;
cf5b0aa0 1577 //
8b63d99c 1578 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
1579 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
1580 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
1581
0b736a46 1582 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
1583 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
cf5b0aa0 1584 pointCov[0]=prot1.GetCov()[3];//simay^2
1585 pointCov[1]=prot1.GetCov()[4];//sigmayz
1586 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
ca58ed4e 1587 if (!track1->Update(pointPos,pointCov)) break;
1588 npoints1++;
1589 npoints2++;
cf5b0aa0 1590 }
cfe2c39a 1591 if (npoints2<npoints/4.) return 0;
1592 AliTrackerBase::PropagateTrackTo(track0,refX,kMass,5.,kTRUE,kMaxSnp);
1593 AliTrackerBase::PropagateTrackTo(track0,refX,kMass,1.,kTRUE,kMaxSnp);
cf5b0aa0 1594 track1->Rotate(track0->GetAlpha());
cfe2c39a 1595 AliTrackerBase::PropagateTrackTo(track1,track0->GetX(),kMass,5.,kFALSE,kMaxSnp);
cf5b0aa0 1596
cad404e1 1597 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
cf5b0aa0 1598 "point0.="<<&pointArray0<< // points
1599 "point1.="<<&pointArray1<< // distorted points
1600 "trackIn.="<<&track<< // original track
1601 "track0.="<<track0<< // fitted track
1602 "track1.="<<track1<< // fitted distorted track
1603 "\n";
be67055b 1604 new(&trackIn) AliExternalTrackParam(*track0);
cf5b0aa0 1605 delete track0;
1606 return track1;
1607}
1608
1609
ffab0c37 1610
1611
1612
1613TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
1614 //
1615 // create the distortion tree on a mesh with granularity given by step
1616 // return the tree with distortions at given position
1617 // Map is created on the mesh with given step size
1618 //
1619 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
1620 Float_t xyz[3];
1621 for (Double_t x= -250; x<250; x+=step){
1622 for (Double_t y= -250; y<250; y+=step){
1623 Double_t r = TMath::Sqrt(x*x+y*y);
1624 if (r<80) continue;
1625 if (r>250) continue;
1626 for (Double_t z= -250; z<250; z+=step){
1627 Int_t roc=(z>0)?0:18;
1628 xyz[0]=x;
1629 xyz[1]=y;
1630 xyz[2]=z;
1631 Double_t phi = TMath::ATan2(y,x);
1632 DistortPoint(xyz,roc);
1633 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1634 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
1635 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
1636 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
1637 Double_t dx = xyz[0]-x;
1638 Double_t dy = xyz[1]-y;
1639 Double_t dz = xyz[2]-z;
1640 Double_t dr=r1-r;
1641 Double_t drphi=(phi1-phi)*r;
1642 (*pcstream)<<"distortion"<<
1643 "x="<<x<< // original position
1644 "y="<<y<<
1645 "z="<<z<<
1646 "r="<<r<<
1647 "phi="<<phi<<
1648 "x1="<<xyz[0]<< // distorted position
1649 "y1="<<xyz[1]<<
1650 "z1="<<xyz[2]<<
1651 "r1="<<r1<<
1652 "phi1="<<phi1<<
1653 //
1654 "dx="<<dx<< // delta position
1655 "dy="<<dy<<
1656 "dz="<<dz<<
1657 "dr="<<dr<<
1658 "drphi="<<drphi<<
1659 "\n";
1660 }
1661 }
1662 }
1663 delete pcstream;
1664 TFile f(Form("correction%s.root",GetName()));
1665 TTree * tree = (TTree*)f.Get("distortion");
1666 TTree * tree2= tree->CopyTree("1");
1667 tree2->SetName(Form("dist%s",GetName()));
1668 tree2->SetDirectory(0);
1669 delete tree;
1670 return tree2;
1671}
1672
1673
1674
be67055b 1675
46e89793 1676void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
be67055b 1677 //
1678 // Make a fit tree:
1679 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1680 // calculates partial distortions
1681 // Partial distortion is stored in the resulting tree
1682 // Output is storred in the file distortion_<dettype>_<partype>.root
1683 // Partial distortion is stored with the name given by correction name
1684 //
1685 //
1686 // Parameters of function:
1687 // input - input tree
cfe2c39a 1688 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex , 3 - TPC-TOF, 4 - TPCTPC track crossing
be67055b 1689 // ppype - parameter type
1690 // corrArray - array with partial corrections
1691 // step - skipe entries - if 1 all entries processed - it is slow
1692 // debug 0 if debug on also space points dumped - it is slow
c9cbd2f2 1693
b322e06a 1694 const Double_t kMaxSnp = 0.85;
cfe2c39a 1695 const Double_t kcutSnp=0.25;
1696 const Double_t kcutTheta=1.;
1697 const Double_t kRadiusTPC=85;
1698 // AliTPCROC *tpcRoc =AliTPCROC::Instance();
1699 //
b322e06a 1700 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1701 // const Double_t kB2C=-0.299792458e-3;
46e89793 1702 const Int_t kMinEntries=20;
cfe2c39a 1703 Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
46e89793 1704 Float_t refX;
1705 Int_t run;
1706 tinput->SetBranchAddress("run",&run);
be67055b 1707 tinput->SetBranchAddress("theta",&theta);
1708 tinput->SetBranchAddress("phi", &phi);
1709 tinput->SetBranchAddress("snp",&snp);
1710 tinput->SetBranchAddress("mean",&mean);
1711 tinput->SetBranchAddress("rms",&rms);
1712 tinput->SetBranchAddress("entries",&entries);
cfe2c39a 1713 tinput->SetBranchAddress("sector",&sector);
1714 tinput->SetBranchAddress("dsec",&dsec);
1715 tinput->SetBranchAddress("refX",&refX);
46e89793 1716 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
be67055b 1717 //
1718 Int_t nentries=tinput->GetEntries();
1719 Int_t ncorr=corrArray->GetEntries();
7f4cb119 1720 Double_t corrections[100]={0}; //
be67055b 1721 Double_t tPar[5];
1722 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
be67055b 1723 Int_t dir=0;
cfe2c39a 1724 if (dtype==5 || dtype==6) dtype=4;
1725 if (dtype==0) { dir=-1;}
1726 if (dtype==1) { dir=1;}
1727 if (dtype==2) { dir=-1;}
1728 if (dtype==3) { dir=1;}
1729 if (dtype==4) { dir=-1;}
be67055b 1730 //
46e89793 1731 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
be67055b 1732 tinput->GetEntry(ientry);
7f4cb119 1733 if (TMath::Abs(snp)>kMaxSnp) continue;
be67055b 1734 tPar[0]=0;
1735 tPar[1]=theta*refX;
cfe2c39a 1736 if (dtype==2) tPar[1]=theta*kRadiusTPC;
be67055b 1737 tPar[2]=snp;
1738 tPar[3]=theta;
4486a91f 1739 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
cfe2c39a 1740 if (dtype==4){
1741 // tracks crossing CE
1742 tPar[1]=0; // track at the CE
1743 //if (TMath::Abs(theta) <0.05) continue; // deep cross
1744 }
1745
1746 if (TMath::Abs(snp) >kcutSnp) continue;
1747 if (TMath::Abs(theta) >kcutTheta) continue;
1748 printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
8b63d99c 1749 Double_t bz=AliTrackerBase::GetBz();
cfe2c39a 1750 if (dtype !=4) { //exclude TPC - for TPC mainly non primary tracks
1751 if (dtype!=2 && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
1752
1753 if (dtype==2 && TMath::Abs(bz)>0.1 ) {
1754 tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);//
1755 // snp at the TPC inner radius in case the vertex match used
1756 }
1757 }
1758 //
4486a91f 1759 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
7f4cb119 1760 AliExternalTrackParam track(refX,phi,tPar,cov);
1761 Double_t xyz[3];
1762 track.GetXYZ(xyz);
1763 Int_t id=0;
46e89793 1764 Double_t pt=1./tPar[4];
7f4cb119 1765 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
cfe2c39a 1766 //if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used
46e89793 1767 Double_t refXD=refX;
be67055b 1768 (*pcstream)<<"fit"<<
46e89793 1769 "run="<<run<< // run number
8b63d99c 1770 "bz="<<bz<< // magnetic filed used
be67055b 1771 "dtype="<<dtype<< // detector match type
1772 "ptype="<<ptype<< // parameter type
1773 "theta="<<theta<< // theta
1774 "phi="<<phi<< // phi
1775 "snp="<<snp<< // snp
1776 "mean="<<mean<< // mean dist value
1777 "rms="<<rms<< // rms
cfe2c39a 1778 "sector="<<sector<<
1779 "dsec="<<dsec<<
46e89793 1780 "refX="<<refXD<< // referece X as double
7f4cb119 1781 "gx="<<xyz[0]<< // global position at reference
1782 "gy="<<xyz[1]<< // global position at reference
1783 "gz="<<xyz[2]<< // global position at reference
1784 "dRrec="<<dRrec<< // delta Radius in reconstruction
46e89793 1785 "pt="<<pt<< // pt
7f4cb119 1786 "id="<<id<< // track id
be67055b 1787 "entries="<<entries;// number of entries in bin
1788 //
cfe2c39a 1789 Bool_t isOK=kTRUE;
46e89793 1790 if (entries<kMinEntries) isOK=kFALSE;
1791 //
cfe2c39a 1792 if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
be67055b 1793 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1794 corrections[icorr]=0;
1795 if (entries>kMinEntries){
1796 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1797 AliExternalTrackParam *trackOut = 0;
1798 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1799 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
cfe2c39a 1800 if (dtype==0) {dir= -1;}
1801 if (dtype==1) {dir= 1;}
1802 if (dtype==2) {dir= -1;}
1803 if (dtype==3) {dir= 1;}
b1f0a2a5 1804 //
7f4cb119 1805 if (trackOut){
cfe2c39a 1806 if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1807 if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE;
1808 if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1809 // trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1810 //
7f4cb119 1811 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1812 delete trackOut;
1813 }else{
1814 corrections[icorr]=0;
cfe2c39a 1815 isOK=kFALSE;
7f4cb119 1816 }
cfe2c39a 1817 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out
be67055b 1818 }
1819 (*pcstream)<<"fit"<<
46e89793 1820 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
be67055b 1821 }
cfe2c39a 1822
1823 if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1824 //
1825 // special case of the TPC tracks crossing the CE
1826 //
1827 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1828 corrections[icorr]=0;
1829 if (entries>kMinEntries){
46e89793 1830 AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
1831 AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet
cfe2c39a 1832 AliExternalTrackParam *trackOut0 = 0;
1833 AliExternalTrackParam *trackOut1 = 0;
1834 //
1835 if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1836 if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1837 if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1838 if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1839 //
1840 if (trackOut0 && trackOut1){
1841 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1842 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1843 if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1844 if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1845 //
1846 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1847 if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1848 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1849 if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE;
1850 if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1851 //
1852 corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1853 corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
46e89793 1854 if (isOK)
1855 if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
1856 (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
1857 (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
1858 (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
1859 (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
1860 (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
1861 ){
1862 isOK=kFALSE;
1863 }
cfe2c39a 1864 delete trackOut0;
46e89793 1865 delete trackOut1;
cfe2c39a 1866 }else{
1867 corrections[icorr]=0;
1868 isOK=kFALSE;
1869 }
1870 //
1871 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out no in lookup
1872 }
cfe2c39a 1873 (*pcstream)<<"fit"<<
46e89793 1874 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
cfe2c39a 1875 }
1876 //
1877 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
be67055b 1878 }
cfe2c39a 1879
1880
be67055b 1881 delete pcstream;
1882}
1883
1884
1885
46e89793 1886void AliTPCCorrection::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1887 //
1888 // Make a fit tree:
1889 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1890 // calculates partial distortions
1891 // Partial distortion is stored in the resulting tree
1892 // Output is storred in the file distortion_<dettype>_<partype>.root
1893 // Partial distortion is stored with the name given by correction name
1894 //
1895 //
1896 // Parameters of function:
1897 // input - input tree
1898 // dtype - distortion type 10 - IROC-OROC
1899 // ppype - parameter type
1900 // corrArray - array with partial corrections
1901 // step - skipe entries - if 1 all entries processed - it is slow
1902 // debug 0 if debug on also space points dumped - it is slow
1903
1904 const Double_t kMaxSnp = 0.8;
1905 const Int_t kMinEntries=200;
1906 // AliTPCROC *tpcRoc =AliTPCROC::Instance();
1907 //
1908 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1909 // const Double_t kB2C=-0.299792458e-3;
1910 Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
1911 Int_t isec1, isec0;
1912 Double_t refXD;
1913 Float_t refX;
1914 Int_t run;
1915 tinput->SetBranchAddress("run",&run);
1916 tinput->SetBranchAddress("theta",&theta);
1917 tinput->SetBranchAddress("phi", &phi);
1918 tinput->SetBranchAddress("snp",&snp);
1919 tinput->SetBranchAddress("mean",&mean);
1920 tinput->SetBranchAddress("rms",&rms);
1921 tinput->SetBranchAddress("entries",&entries);
1922 tinput->SetBranchAddress("sector",&sector);
1923 tinput->SetBranchAddress("dsec",&dsec);
1924 tinput->SetBranchAddress("refX",&refXD);
1925 tinput->SetBranchAddress("z",&globalZ);
1926 tinput->SetBranchAddress("isec0",&isec0);
1927 tinput->SetBranchAddress("isec1",&isec1);
1928 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
1929 //
1930 Int_t nentries=tinput->GetEntries();
1931 Int_t ncorr=corrArray->GetEntries();
1932 Double_t corrections[100]={0}; //
1933 Double_t tPar[5];
1934 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1935 Int_t dir=0;
1936 //
1937 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1938 tinput->GetEntry(ientry);
1939 refX=refXD;
1940 Int_t id=-1;
1941 if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1; // IROC-OROC - opposite side
1942 if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2; // IROC-OROC - same side
1943 if (dtype==10 && id==-1) continue;
1944 //
1945 dir=-1;
1946 tPar[0]=0;
1947 tPar[1]=globalZ;
1948 tPar[2]=snp;
1949 tPar[3]=theta;
1950 tPar[4]=(gRandom->Rndm()-0.1)*0.2; //
1951 Double_t pt=1./tPar[4];
1952 //
1953 printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1954 Double_t bz=AliTrackerBase::GetBz();
1955 AliExternalTrackParam track(refX,phi,tPar,cov);
1956 Double_t xyz[3],xyzIn[3],xyzOut[3];
1957 track.GetXYZ(xyz);
1958 track.GetXYZAt(85,bz,xyzIn);
1959 track.GetXYZAt(245,bz,xyzOut);
1960 Double_t phiIn = TMath::ATan2(xyzIn[1],xyzIn[0]);
1961 Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
1962 Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
1963 Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
1964 Int_t sectorIn = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
1965 Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
1966 //
1967 Bool_t isOK=kTRUE;
1968 if (sectorIn!=sectorOut) isOK=kFALSE; // requironment - cluster in the same sector
1969 if (sectorIn!=sectorRef) isOK=kFALSE; // requironment - cluster in the same sector
1970 if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE; // requironment - minimal amount of tracks in bin
1971 // Do downscale
1972 if (TMath::Abs(theta)>1) isOK=kFALSE;
1973 //
1974 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1975 //
1976 (*pcstream)<<"fit"<<
1977 "run="<<run<< //run
1978 "bz="<<bz<< // magnetic filed used
1979 "dtype="<<dtype<< // detector match type
1980 "ptype="<<ptype<< // parameter type
1981 "theta="<<theta<< // theta
1982 "phi="<<phi<< // phi
1983 "snp="<<snp<< // snp
1984 "mean="<<mean<< // mean dist value
1985 "rms="<<rms<< // rms
1986 "sector="<<sector<<
1987 "dsec="<<dsec<<
1988 "refX="<<refXD<< // referece X
1989 "gx="<<xyz[0]<< // global position at reference
1990 "gy="<<xyz[1]<< // global position at reference
1991 "gz="<<xyz[2]<< // global position at reference
1992 "dRrec="<<dRrec<< // delta Radius in reconstruction
1993 "pt="<<pt<< //pt
1994 "id="<<id<< // track id
1995 "entries="<<entries;// number of entries in bin
1996 //
1997 AliExternalTrackParam *trackOut0 = 0;
1998 AliExternalTrackParam *trackOut1 = 0;
1999 AliExternalTrackParam *ptrackIn0 = 0;
2000 AliExternalTrackParam *ptrackIn1 = 0;
2001
2002 for (Int_t icorr=0; icorr<ncorr; icorr++) {
2003 //
2004 // special case of the TPC tracks crossing the CE
2005 //
2006 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
2007 corrections[icorr]=0;
2008 if (entries>kMinEntries &&isOK){
2009 AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
2010 AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
2011 ptrackIn1=&trackIn0;
2012 ptrackIn0=&trackIn1;
2013 //
2014 if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
2015 if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
2016 if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
2017 if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
2018 //
2019 if (trackOut0 && trackOut1){
2020 //
2021 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp)) isOK=kFALSE;
2022 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2023 // rotate all tracks to the same frame
2024 if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2025 if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2026 if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2027 //
2028 if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2029 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2030 if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2031 //
2032 corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
2033 corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
2034 (*pcstream)<<"fitDebug"<< // just to debug the correction
2035 "mean="<<mean<<
2036 "pIn0.="<<ptrackIn0<<
2037 "pIn1.="<<ptrackIn1<<
2038 "pOut0.="<<trackOut0<<
2039 "pOut1.="<<trackOut1<<
2040 "refX="<<refXD<<
2041 "\n";
2042 delete trackOut0;
2043 delete trackOut1;
2044 }else{
2045 corrections[icorr]=0;
2046 isOK=kFALSE;
2047 }
2048 }
2049 (*pcstream)<<"fit"<<
2050 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
2051 }
2052 //
2053 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
2054 }
2055 delete pcstream;
2056}
2057
2058
2059
2060void AliTPCCorrection::MakeLaserDistortionTreeOld(TTree* tree, TObjArray *corrArray, Int_t itype){
7f4cb119 2061 //
2062 // Make a laser fit tree for global minimization
2063 //
2064 const Double_t cutErrY=0.1;
2065 const Double_t cutErrZ=0.1;
2066 const Double_t kEpsilon=0.00000001;
46e89793 2067 const Double_t kMaxDist=1.; // max distance - space correction
2068 const Double_t kMaxRMS=0.05; // max distance -between point and local mean
7f4cb119 2069 TVectorD *vecdY=0;
2070 TVectorD *vecdZ=0;
2071 TVectorD *veceY=0;
2072 TVectorD *veceZ=0;
2073 AliTPCLaserTrack *ltr=0;
2074 AliTPCLaserTrack::LoadTracks();
2075 tree->SetBranchAddress("dY.",&vecdY);
2076 tree->SetBranchAddress("dZ.",&vecdZ);
2077 tree->SetBranchAddress("eY.",&veceY);
2078 tree->SetBranchAddress("eZ.",&veceZ);
2079 tree->SetBranchAddress("LTr.",&ltr);
2080 Int_t entries= tree->GetEntries();
cfe2c39a 2081 TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
7f4cb119 2082 Double_t bz=AliTrackerBase::GetBz();
2083 //
2084
2085 for (Int_t ientry=0; ientry<entries; ientry++){
2086 tree->GetEntry(ientry);
2087 if (!ltr->GetVecGX()){
2088 ltr->UpdatePoints();
2089 }
2090 TVectorD * delta= (itype==0)? vecdY:vecdZ;
2091 TVectorD * err= (itype==0)? veceY:veceZ;
46e89793 2092 TLinearFitter fitter(2,"pol1");
2093 for (Int_t iter=0; iter<2; iter++){
2094 Double_t kfit0=0, kfit1=0;
2095 Int_t npoints=fitter.GetNpoints();
2096 if (npoints>80){
2097 fitter.Eval();
2098 kfit0=fitter.GetParameter(0);
2099 kfit1=fitter.GetParameter(1);
2100 }
2101 for (Int_t irow=0; irow<159; irow++){
2102 Bool_t isOK=kTRUE;
2103 Int_t isOKF=0;
2104 Int_t nentries = 1000;
2105 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
2106 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
2107 Int_t dtype=5;
2108 Double_t array[10];
2109 Int_t first3=TMath::Max(irow-3,0);
2110 Int_t last3 =TMath::Min(irow+3,159);
2111 Int_t counter=0;
2112 if ((*ltr->GetVecSec())[irow]>=0 && err) {
2113 for (Int_t jrow=first3; jrow<=last3; jrow++){
2114 if ((*ltr->GetVecSec())[irow]!= (*ltr->GetVecSec())[jrow]) continue;
2115 if ((*err)[jrow]<kEpsilon) continue;
2116 array[counter]=(*delta)[jrow];
2117 counter++;
2118 }
2119 }
2120 Double_t rms3 = 0;
2121 Double_t mean3 = 0;
2122 if (counter>2){
2123 rms3 = TMath::RMS(counter,array);
2124 mean3 = TMath::Mean(counter,array);
2125 }else{
2126 isOK=kFALSE;
2127 }
2128 Double_t phi =(*ltr->GetVecPhi())[irow];
2129 Double_t theta =ltr->GetTgl();
2130 Double_t mean=delta->GetMatrixArray()[irow];
2131 Double_t gx=0,gy=0,gz=0;
2132 Double_t snp = (*ltr->GetVecP2())[irow];
2133 Double_t dRrec=0;
2134 // Double_t rms = err->GetMatrixArray()[irow];
cfe2c39a 2135 //
46e89793 2136 gx = (*ltr->GetVecGX())[irow];
2137 gy = (*ltr->GetVecGY())[irow];
2138 gz = (*ltr->GetVecGZ())[irow];
2139 //
2140 // get delta R used in reconstruction
2141 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
2142 AliTPCCorrection * correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
2143 // const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
2144 //Double_t xyz0[3]={gx,gy,gz};
2145 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
2146 Double_t fphi = TMath::ATan2(gy,gx);
2147 Double_t fsector = 9.*fphi/TMath::Pi();
2148 if (fsector<0) fsector+=18;
2149 Double_t dsec = fsector-Int_t(fsector)-0.5;
2150 Double_t refX=0;
2151 Int_t id= ltr->GetId();
2152 Double_t pt=0;
2153 //
2154 if (1 && oldR>1) {
2155 Float_t xyz1[3]={gx,gy,gz};
2156 Int_t sector=(gz>0)?0:18;
2157 correction->CorrectPoint(xyz1, sector);
2158 refX=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
2159 dRrec=oldR-refX;
2160 }
2161 if (TMath::Abs(rms3)>kMaxRMS) isOK=kFALSE;
2162 if (TMath::Abs(mean-mean3)>kMaxRMS) isOK=kFALSE;
2163 if (counter<4) isOK=kFALSE;
2164 if (npoints<90) isOK=kFALSE;
2165 if (isOK){
2166 fitter.AddPoint(&refX,mean);
7f4cb119 2167 }
46e89793 2168 Double_t deltaF=kfit0+kfit1*refX;
2169 if (iter==1){
2170 (*pcstream)<<"fitFull"<< // dumpe also intermediate results
2171 "bz="<<bz<< // magnetic filed used
2172 "dtype="<<dtype<< // detector match type
2173 "ptype="<<itype<< // parameter type
2174 "theta="<<theta<< // theta
2175 "phi="<<phi<< // phi
2176 "snp="<<snp<< // snp
2177 "mean="<<mean3<< // mean dist value
2178 "rms="<<rms3<< // rms
2179 "deltaF="<<deltaF<<
2180 "npoints="<<npoints<< //number of points
2181 "mean3="<<mean3<< // mean dist value
2182 "rms3="<<rms3<< // rms
2183 "counter="<<counter<<
2184 "sector="<<fsector<<
2185 "dsec="<<dsec<<
2186 //
2187 "refX="<<refX<< // reference radius
2188 "gx="<<gx<< // global position
2189 "gy="<<gy<< // global position
2190 "gz="<<gz<< // global position
2191 "dRrec="<<dRrec<< // delta Radius in reconstruction
2192 "id="<<id<< //bundle
2193 "entries="<<nentries<<// number of entries in bin
2194 "\n";
2195 }
2196 if (iter==1) (*pcstream)<<"fit"<< // dump valus for fit
2197 "bz="<<bz<< // magnetic filed used
2198 "dtype="<<dtype<< // detector match type
2199 "ptype="<<itype<< // parameter type
2200 "theta="<<theta<< // theta
2201 "phi="<<phi<< // phi
2202 "snp="<<snp<< // snp
2203 "mean="<<mean3<< // mean dist value
2204 "rms="<<rms3<< // rms
2205 "sector="<<fsector<<
2206 "dsec="<<dsec<<
2207 //
2208 "refX="<<refX<< // reference radius
2209 "gx="<<gx<< // global position
2210 "gy="<<gy<< // global position
2211 "gz="<<gz<< // global position
2212 "dRrec="<<dRrec<< // delta Radius in reconstruction
2213 "pt="<<pt<< //pt
2214 "id="<<id<< //bundle
2215 "entries="<<nentries;// number of entries in bin
2216 //
2217 //
2218 Double_t ky = TMath::Tan(TMath::ASin(snp));
2219 Int_t ncorr = corrArray->GetEntries();
2220 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
2221 Double_t phi0 = TMath::ATan2(gy,gx);
2222 Double_t distortions[1000]={0};
2223 Double_t distortionsR[1000]={0};
2224 if (iter==1){
2225 for (Int_t icorr=0; icorr<ncorr; icorr++) {
2226 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
2227 Float_t distPoint[3]={gx,gy,gz};
2228 Int_t sector= (gz>0)? 0:18;
2229 if (r0>80){
2230 corr->DistortPoint(distPoint, sector);
2231 }
2232 // Double_t value=distPoint[2]-gz;
2233 if (itype==0 && r0>1){
2234 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
2235 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
2236 Double_t drphi= r0*(phi1-phi0);
2237 Double_t dr = r1-r0;
2238 distortions[icorr] = drphi-ky*dr;
2239 distortionsR[icorr] = dr;
2240 }
2241 if (TMath::Abs(distortions[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE; }
2242 if (TMath::Abs(distortionsR[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE;}
2243 (*pcstream)<<"fit"<<
2244 Form("%s=",corr->GetName())<<distortions[icorr]; // dump correction value
2245 }
2246 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
7f4cb119 2247 }
7f4cb119 2248 }
7f4cb119 2249 }
2250 }
2251 delete pcstream;
2252}
2253
2254
be67055b 2255
97d17739 2256void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type, Int_t integ){
cfe2c39a 2257 //
2258 // make a distortion map out ou fthe residual histogram
2259 // Results are written to the debug streamer - pcstream
2260 // Parameters:
2261 // his0 - input (4D) residual histogram
2262 // pcstream - file to write the tree
2263 // run - run number
2264 // refX - track matching reference X
2265 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
2266 // THnSparse axes:
2267 // OBJ: TAxis #Delta #Delta
2268 // OBJ: TAxis tanTheta tan(#Theta)
2269 // OBJ: TAxis phi #phi
2270 // OBJ: TAxis snp snp
2271
2272 // marian.ivanov@cern.ch
2273 const Int_t kMinEntries=10;
2274 Double_t bz=AliTrackerBase::GetBz();
2275 Int_t idim[4]={0,1,2,3};
2276 //
2277 //
2278 //
2279 Int_t nbins3=his0->GetAxis(3)->GetNbins();
2280 Int_t first3=his0->GetAxis(3)->GetFirst();
2281 Int_t last3 =his0->GetAxis(3)->GetLast();
2282 //
2283 for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - local angle
97d17739 2284 his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3));
cfe2c39a 2285 Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3);
2286 THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3
2287 //
2288 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
2289 Int_t first2 = his3->GetAxis(2)->GetFirst();
2290 Int_t last2 = his3->GetAxis(2)->GetLast();
2291 //
2292 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi
97d17739 2293 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2));
cfe2c39a 2294 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
2295 THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
2296 Int_t nbins1 = his2->GetAxis(1)->GetNbins();
2297 Int_t first1 = his2->GetAxis(1)->GetFirst();
2298 Int_t last1 = his2->GetAxis(1)->GetLast();
2299 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
2300 //
2301 Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1);
2302 his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
2303 if (TMath::Abs(x1)<0.1){
2304 if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1));
2305 if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1));
2306 }
2307 if (TMath::Abs(x1)<0.06){
2308 his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
2309 }
2310 TH1 * hisDelta = his2->Projection(0);
2311 //
2312 Double_t entries = hisDelta->GetEntries();
2313 Double_t mean=0, rms=0;
2314 if (entries>kMinEntries){
2315 mean = hisDelta->GetMean();
2316 rms = hisDelta->GetRMS();
2317 }
2318 Double_t sector = 9.*x2/TMath::Pi();
2319 if (sector<0) sector+=18;
2320 Double_t dsec = sector-Int_t(sector)-0.5;
2321 Double_t z=refX*x1;
2322 (*pcstream)<<hname<<
2323 "run="<<run<<
2324 "bz="<<bz<<
2325 "theta="<<x1<<
2326 "phi="<<x2<<
2327 "z="<<z<< // dummy z
2328 "snp="<<x3<<
2329 "entries="<<entries<<
2330 "mean="<<mean<<
2331 "rms="<<rms<<
2332 "refX="<<refX<< // track matching refernce plane
2333 "type="<<type<< //
2334 "sector="<<sector<<
2335 "dsec="<<dsec<<
2336 "\n";
2337 delete hisDelta;
02cd5ade 2338 //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
cfe2c39a 2339 }
2340 delete his2;
2341 }
2342 delete his3;
2343 }
2344}
2345
2346
2347
2348
2349void AliTPCCorrection::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){
8b63d99c 2350 //
2351 // make a distortion map out ou fthe residual histogram
2352 // Results are written to the debug streamer - pcstream
2353 // Parameters:
2354 // his0 - input (4D) residual histogram
2355 // pcstream - file to write the tree
2356 // run - run number
cfe2c39a 2357 // refX - track matching reference X
2358 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
8b63d99c 2359 // marian.ivanov@cern.ch
cfe2c39a 2360 //
2361 // Histo axeses
2362 // Collection name='TObjArray', class='TObjArray', size=16
2363 // 0. OBJ: TAxis #Delta #Delta
2364 // 1. OBJ: TAxis N_{cl} N_{cl}
2365 // 2. OBJ: TAxis dca_{r} (cm) dca_{r} (cm)
2366 // 3. OBJ: TAxis z (cm) z (cm)
2367 // 4. OBJ: TAxis sin(#phi) sin(#phi)
2368 // 5. OBJ: TAxis tan(#theta) tan(#theta)
2369 // 6. OBJ: TAxis 1/pt (1/GeV) 1/pt (1/GeV)
2370 // 7. OBJ: TAxis pt (GeV) pt (GeV)
2371 // 8. OBJ: TAxis alpha alpha
2372 const Int_t kMinEntries=10;
2373 //
2374 // 1. make default selections
2375 //
2376 TH1 * hisDelta=0;
2377 Int_t idim0[4]={0 , 5, 8, 3}; // delta, theta, alpha, z
2378 hisInput->GetAxis(1)->SetRangeUser(110,190); //long tracks
2379 hisInput->GetAxis(2)->SetRangeUser(-10,35); //tracks close to beam pipe
2380 hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance
2381 hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks"
2382 hisDelta= hisInput->Projection(0);
2383 hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS());
2384 delete hisDelta;
2385 THnSparse *his0= hisInput->Projection(4,idim0);
2386 //
2387 // 2. Get mean in diferent bins
2388 //
8b63d99c 2389 Int_t nbins1=his0->GetAxis(1)->GetNbins();
2390 Int_t first1=his0->GetAxis(1)->GetFirst();
2391 Int_t last1 =his0->GetAxis(1)->GetLast();
2392 //
2393 Double_t bz=AliTrackerBase::GetBz();
cfe2c39a 2394 Int_t idim[4]={0,1, 2, 3}; // delta, theta,alpha,z
2395 //
2396 for (Int_t ibin1=first1; ibin1<=last1; ibin1++){ //axis 1 - theta
2397 //
2398 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
2399 his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
8b63d99c 2400 //
8b63d99c 2401 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
2402 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
2403 Int_t first3 = his1->GetAxis(3)->GetFirst();
2404 Int_t last3 = his1->GetAxis(3)->GetLast();
2405 //
cfe2c39a 2406 for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - z at "vertex"
8b63d99c 2407 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
2408 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
2409 if (ibin3<first3) {
2410 his1->GetAxis(3)->SetRangeUser(-1,1);
2411 x3=0;
2412 }
2413 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
2414 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
2415 Int_t first2 = his3->GetAxis(2)->GetFirst();
2416 Int_t last2 = his3->GetAxis(2)->GetLast();
2417 //
cfe2c39a 2418 for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){
8b63d99c 2419 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
2420 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
cfe2c39a 2421 hisDelta = his3->Projection(0);
8b63d99c 2422 //
2423 Double_t entries = hisDelta->GetEntries();
2424 Double_t mean=0, rms=0;
2425 if (entries>kMinEntries){
2426 mean = hisDelta->GetMean();
2427 rms = hisDelta->GetRMS();
2428 }
cfe2c39a 2429 Double_t sector = 9.*x2/TMath::Pi();
2430 if (sector<0) sector+=18;
2431 Double_t dsec = sector-Int_t(sector)-0.5;
2432 Double_t snp=0; // dummy snp - equal 0
8b63d99c 2433 (*pcstream)<<hname<<
2434 "run="<<run<<
cfe2c39a 2435 "bz="<<bz<< // magnetic field
2436 "theta="<<x1<< // theta
2437 "phi="<<x2<< // phi (alpha)
2438 "z="<<x3<< // z at "vertex"
2439 "snp="<<snp<< // dummy snp
2440 "entries="<<entries<< // entries in bin
2441 "mean="<<mean<< // mean
8b63d99c 2442 "rms="<<rms<<
cfe2c39a 2443 "refX="<<refX<< // track matching refernce plane
2444 "type="<<type<< // parameter type
2445 "sector="<<sector<< // sector
2446 "dsec="<<dsec<< // dummy delta sector
8b63d99c 2447 "\n";
2448 delete hisDelta;
2449 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
2450 }
2451 delete his3;
2452 }
2453 delete his1;
2454 }
cfe2c39a 2455 delete his0;
2456}
2457
2458
2459
2460void AliTPCCorrection::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type){
2461 //
2462 // make a distortion map out of the residual histogram
2463 // Results are written to the debug streamer - pcstream
2464 // Parameters:
2465 // his0 - input (4D) residual histogram
2466 // pcstream - file to write the tree
2467 // run - run number
2468 // type - 0- y 1-z,2 -snp, 3-theta
2469 // marian.ivanov@cern.ch
2470
2471 //Collection name='TObjArray', class='TObjArray', size=16
2472 //0 OBJ: TAxis delta delta
2473 //1 OBJ: TAxis phi phi
2474 //2 OBJ: TAxis localX localX
2475 //3 OBJ: TAxis kY kY
2476 //4 OBJ: TAxis kZ kZ
2477 //5 OBJ: TAxis is1 is1
2478 //6 OBJ: TAxis is0 is0
2479 //7. OBJ: TAxis z z
2480 //8. OBJ: TAxis IsPrimary IsPrimary
2481
2482 const Int_t kMinEntries=10;
2483 THnSparse * hisSector0=0;
2484 TH1 * htemp=0; // histogram to calculate mean value of parameter
2485 Double_t bz=AliTrackerBase::GetBz();
2486
2487 //
2488 // Loop over pair of sector:
2489 // isPrim - 8 ==> 8
2490 // isec0 - 6 ==> 7
2491 // isec1 - 5 ==> 6
2492 // refX - 2 ==> 5
2493 //
2494 // phi - 1 ==> 4
2495 // z - 7 ==> 3
2496 // snp - 3 ==> 2
2497 // theta- 4 ==> 1
2498 // 0 ==> 0;
2499 for (Int_t isec0=0; isec0<72; isec0++){
2500 Int_t index0[9]={0, 4, 3, 7, 1, 2, 5, 6,8}; //regroup indeces
2501 //
2502 //hisInput->GetAxis(8)->SetRangeUser(-0.1,0.4); // select secondaries only ? - get out later ?
2503 hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1);
2504 hisSector0=hisInput->Projection(7,index0);
2505 //
2506 //
2507 for (Int_t isec1=isec0+1; isec1<72; isec1++){
2508 //if (isec1!=isec0+36) continue;
2509 if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue;
2510 printf("Sectors %d\t%d\n",isec1,isec0);
2511 hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1);
2512 TH1 * hisX=hisSector0->Projection(5);
2513 Double_t refX= hisX->GetMean();
2514 delete hisX;
2515 TH1 *hisDelta=hisSector0->Projection(0);
2516 Double_t dmean = hisDelta->GetMean();
2517 Double_t drms = hisDelta->GetRMS();
2518 hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms);
2519 delete hisDelta;
2520 //
2521 // 1. make default selections
2522 //
2523 Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi }
2524 THnSparse *hisSector1= hisSector0->Projection(5,idim0);
2525 //
2526 // 2. Get mean in diferent bins
2527 //
2528 Int_t idim[5]={0, 1, 2, 3, 4}; // {delta, theta-1,snp-2 ,z-3, phi-4}
2529 //
2530 // Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins();
2531 Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst();
2532 Int_t lastPhi =hisSector1->GetAxis(4)->GetLast();
2533 //
2534 for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=1){ //axis 4 - phi
2535 //
2536 // Phi loop
2537 //
2538 Double_t xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi);
2539 Double_t psec = (9*xPhi/TMath::Pi());
2540 if (psec<0) psec+=18;
2541 Bool_t isOK0=kFALSE;
2542 Bool_t isOK1=kFALSE;
2543 if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.) isOK0=kTRUE;
2544 if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.) isOK1=kTRUE;
2545 if (!isOK0) continue;
2546 if (!isOK1) continue;
2547 //
2548 hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi));
2549 if (isec1!=isec0+36) {
2550 hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi));
2551 }
2552 //
2553 htemp = hisSector1->Projection(4);
2554 xPhi=htemp->GetMean();
2555 delete htemp;
2556 THnSparse * hisPhi = hisSector1->Projection(4,idim);
2557 //Int_t nbinsZ = hisPhi->GetAxis(3)->GetNbins();
2558 Int_t firstZ = hisPhi->GetAxis(3)->GetFirst();
2559 Int_t lastZ = hisPhi->GetAxis(3)->GetLast();
2560 //
2561 for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=1){ // axis 3 - z
2562 //
2563 // Z loop
2564 //
2565 hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ));
2566 if (isec1!=isec0+36) {
2567 hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ));
2568 }
2569 htemp = hisPhi->Projection(3);
2570 Double_t xZ= htemp->GetMean();
2571 delete htemp;
2572 THnSparse * hisZ= hisPhi->Projection(3,idim);
2573 //projected histogram according selection 3 -z
2574 //
2575 //
2576 //Int_t nbinsSnp = hisZ->GetAxis(2)->GetNbins();
2577 Int_t firstSnp = hisZ->GetAxis(2)->GetFirst();
2578 Int_t lastSnp = hisZ->GetAxis(2)->GetLast();
2579 for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){ // axis 2 - snp
2580 //
2581 // Snp loop
2582 //
2583 hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp));
2584 if (isec1!=isec0+36) {
2585 hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp));
2586 }
2587 htemp = hisZ->Projection(2);
2588 Double_t xSnp= htemp->GetMean();
2589 delete htemp;
2590 THnSparse * hisSnp= hisZ->Projection(2,idim);
2591 //projected histogram according selection 2 - snp
2592
2593 //Int_t nbinsTheta = hisSnp->GetAxis(1)->GetNbins();
2594 Int_t firstTheta = hisSnp->GetAxis(1)->GetFirst();
2595 Int_t lastTheta = hisSnp->GetAxis(1)->GetLast();
2596 //
2597 for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){ // axis1 theta
2598
2599
2600 hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta));
2601 if (isec1!=isec0+36) {
2602 hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta));
2603 }
2604 htemp = hisSnp->Projection(1);
2605 Double_t xTheta=htemp->GetMean();
2606 delete htemp;
2607 hisDelta = hisSnp->Projection(0);
2608 //
2609 Double_t entries = hisDelta->GetEntries();
2610 Double_t mean=0, rms=0;
2611 if (entries>kMinEntries){
2612 mean = hisDelta->GetMean();
2613 rms = hisDelta->GetRMS();
2614 }
2615 Double_t sector = 9.*xPhi/TMath::Pi();
2616 if (sector<0) sector+=18;
2617 Double_t dsec = sector-Int_t(sector)-0.5;
2618 Int_t dtype=1; // TPC alignment type
2619 (*pcstream)<<hname<<
2620 "run="<<run<<
2621 "bz="<<bz<< // magnetic field
2622 "ptype="<<type<< // parameter type
2623 "dtype="<<dtype<< // parameter type
2624 "isec0="<<isec0<< // sector 0
2625 "isec1="<<isec1<< // sector 1
2626 "sector="<<sector<< // sector as float
2627 "dsec="<<dsec<< // delta sector
2628 //
2629 "theta="<<xTheta<< // theta
2630 "phi="<<xPhi<< // phi (alpha)
2631 "z="<<xZ<< // z
2632 "snp="<<xSnp<< // snp
2633 //
2634 "entries="<<entries<< // entries in bin
2635 "mean="<<mean<< // mean
2636 "rms="<<rms<< // rms
2637 "refX="<<refX<< // track matching reference plane
2638 "\n";
2639 delete hisDelta;
2640 printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\n",isec0, isec1, xPhi,xZ,xSnp, xTheta, entries,mean);
2641 //
2642 }//ibinTheta
2643 delete hisSnp;
2644 } //ibinSnp
2645 delete hisZ;
2646 }//ibinZ
2647 delete hisPhi;
2648 }//ibinPhi
2649 delete hisSector1;
2650 }//isec1
2651 delete hisSector0;
2652 }//isec0
8b63d99c 2653}
2654
2655
2656
2657
2658
cfe2c39a 2659
2660
ffab0c37 2661void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
2662 //
2663 // Store object in the OCDB
2664 // By default the object is stored in the current directory
2665 // default comment consit of user name and the date
2666 //
2667 TString ocdbStorage="";
2668 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
2669 AliCDBMetaData *metaData= new AliCDBMetaData();
2670 metaData->SetObjectClassName("AliTPCCorrection");
2671 metaData->SetResponsible("Marian Ivanov");
2672 metaData->SetBeamPeriod(1);
2673 metaData->SetAliRootVersion("05-25-01"); //root version
2674 TString userName=gSystem->GetFromPipe("echo $USER");
2675 TString date=gSystem->GetFromPipe("date");
2676
2677 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
2678 if (comment) metaData->SetComment(comment);
2679 AliCDBId* id1=NULL;
2680 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
2681 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
2682 gStorage->Put(this, (*id1), metaData);
2683}
2684
ca58ed4e 2685
7d85e147 2686void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
c9cbd2f2 2687 //
2688 // Fast method to simulate the influence of the given distortion on the vertex reconstruction
2689 //
ca58ed4e 2690
c9cbd2f2 2691 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
2692 if (!magF) AliError("Magneticd field - not initialized");
2693 Double_t bz = magF->SolenoidField(); //field in kGauss
9f3b99e2 2694 printf("bz: %f\n",bz);
c9cbd2f2 2695 AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
ca58ed4e 2696
c9cbd2f2 2697 TObjArray aTrk; // Original Track array of Aside
2698 TObjArray daTrk; // Distorted Track array of A side
2699 UShort_t *aId = new UShort_t[nTracks]; // A side Track ID
2700 TObjArray cTrk;
2701 TObjArray dcTrk;
2702 UShort_t *cId = new UShort_t [nTracks];
2703 Int_t id=0;
ca58ed4e 2704 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
7d85e147 2705 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
ca58ed4e 2706 fpt.SetParameters(7.24,0.120);
2707 fpt.SetNpx(10000);
2708 for(Int_t nt=0; nt<nTracks; nt++){
2709 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
7d85e147 2710 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
c9cbd2f2 2711 Double_t pt = fpt.GetRandom(); // momentum for f1
2712 // printf("phi %lf eta %lf pt %lf\n",phi,eta,pt);
ca58ed4e 2713 Short_t sign=1;
2714 if(gRandom->Rndm() < 0.5){
2715 sign =1;
2716 }else{
2717 sign=-1;
2718 }
2719
2720 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
2721 Double_t pxyz[3];
2722 pxyz[0]=pt*TMath::Cos(phi);
2723 pxyz[1]=pt*TMath::Sin(phi);
2724 pxyz[2]=pt*TMath::Tan(theta);
2725 Double_t cv[21]={0};
2726 AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
2727
2728 Double_t refX=1.;
2729 Int_t dir=-1;
2730 AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir, NULL);
2731 if (!td) continue;
2732 if (pcstream) (*pcstream)<<"track"<<
2733 "eta="<<eta<<
2734 "theta="<<theta<<
2735 "tOrig.="<<t<<
2736 "td.="<<td<<
2737 "\n";
7d85e147 2738 if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
ca58ed4e 2739 if (td){
c9cbd2f2 2740 daTrk.AddLast(td);
2741 aTrk.AddLast(t);
2742 Int_t nn=aTrk.GetEntriesFast();
2743 aId[nn]=id;
ca58ed4e 2744 }
7d85e147 2745 }else if(( eta<-0.07 )&&( eta>-etaCuts )){
ca58ed4e 2746 if (td){
c9cbd2f2 2747 dcTrk.AddLast(td);
2748 cTrk.AddLast(t);
2749 Int_t nn=cTrk.GetEntriesFast();
2750 cId[nn]=id;
ca58ed4e 2751 }
2752 }
c9cbd2f2 2753 id++;
ca58ed4e 2754 }// end of track loop
2755
2756 vertexer->SetTPCMode();
2757 vertexer->SetConstraintOff();
2758
c9cbd2f2 2759 aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));
2760 avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
2761 cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));
2762 cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
ca58ed4e 2763 if (pcstream) (*pcstream)<<"vertex"<<
2764 "x="<<orgVertex[0]<<
2765 "y="<<orgVertex[1]<<
2766 "z="<<orgVertex[2]<<
2767 "av.="<<&aV<< // distorted vertex A side
2768 "cv.="<<&cV<< // distroted vertex C side
2769 "avO.="<<&avOrg<< // original vertex A side
2770 "cvO.="<<&cvOrg<<
2771 "\n";
c9cbd2f2 2772 delete []aId;
2773 delete []cId;
ca58ed4e 2774}
f1817479 2775
2776void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
2777 //
2778 // make correction available for visualization using
2779 // TFormula, TFX and TTree::Draw
2780 // important in order to check corrections and also compute dervied variables
2781 // e.g correction partial derivatives
2782 //
2783 // NOTE - class is not owner of correction
2784 //
cfe2c39a 2785 if (!fgVisualCorrection) fgVisualCorrection=new TObjArray(10000);
2786 if (position>=fgVisualCorrection->GetEntriesFast())
2787 fgVisualCorrection->Expand((position+10)*2);
f1817479 2788 fgVisualCorrection->AddAt(corr, position);
2789}
2790
2791
2792
287fbdfa 2793Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t r, Double_t kZ, Int_t axisType, Int_t corrType){
f1817479 2794 //
2795 // calculate the correction at given position - check the geffCorr
2796 //
cfe2c39a 2797 // corrType return values
2798 // 0 - delta R
2799 // 1 - delta RPhi
2800 // 2 - delta Z
2801 // 3 - delta RPHI
2802 //
f1817479 2803 if (!fgVisualCorrection) return 0;
2804 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
2805 if (!corr) return 0;
25732bff 2806
f1817479 2807 Double_t phi=sector*TMath::Pi()/9.;
287fbdfa 2808 Double_t gx = r*TMath::Cos(phi);
2809 Double_t gy = r*TMath::Sin(phi);
2810 Double_t gz = r*kZ;
f1817479 2811 Int_t nsector=(gz>0) ? 0:18;
2812 //
2813 //
2814 //
2815 Float_t distPoint[3]={gx,gy,gz};
2816 corr->DistortPoint(distPoint, nsector);
2817 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
2818 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
2819 Double_t phi0=TMath::ATan2(gy,gx);
2820 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
2821 if (axisType==0) return r1-r0;
2822 if (axisType==1) return (phi1-phi0)*r0;
2823 if (axisType==2) return distPoint[2]-gz;
cfe2c39a 2824 if (axisType==3) return (TMath::Cos(phi)*(distPoint[0]-gx)+ TMath::Cos(phi)*(distPoint[1]-gy));
f1817479 2825 return phi1-phi0;
2826}
2827
2828Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
2829 //
2830 // return correction at given x,y,z
2831 //
2832 if (!fgVisualCorrection) return 0;
2833 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
2834 if (!corr) return 0;
2835 Double_t phi0= TMath::ATan2(gy,gx);
2836 Int_t nsector=(gz>0) ? 0:18;
2837 Float_t distPoint[3]={gx,gy,gz};
2838 corr->DistortPoint(distPoint, nsector);
2839 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
2840 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
2841 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
2842 if (axisType==0) return r1-r0;
2843 if (axisType==1) return (phi1-phi0)*r0;
2844 if (axisType==2) return distPoint[2]-gz;
2845 return phi1-phi0;
2846}
46e89793 2847
2848
2849
2850
2851
284418bc 2852void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){
46e89793 2853 //
2854 // Make a laser fit tree for global minimization
2855 //
2856 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
2857 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
2858 if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
2859 correction->AddVisualCorrection(correction,0); //register correction
2860
284418bc 2861 // AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2862 //AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
46e89793 2863 //
2864 const Double_t cutErrY=0.05;
2865 const Double_t kSigmaCut=4;
2866 // const Double_t cutErrZ=0.03;
2867 const Double_t kEpsilon=0.00000001;
284418bc 2868 // const Double_t kMaxDist=1.; // max distance - space correction
46e89793 2869 TVectorD *vecdY=0;
2870 TVectorD *vecdZ=0;
2871 TVectorD *veceY=0;
2872 TVectorD *veceZ=0;
2873 AliTPCLaserTrack *ltr=0;
2874 AliTPCLaserTrack::LoadTracks();
2875 tree->SetBranchAddress("dY.",&vecdY);
2876 tree->SetBranchAddress("dZ.",&vecdZ);
2877 tree->SetBranchAddress("eY.",&veceY);
2878 tree->SetBranchAddress("eZ.",&veceZ);
2879 tree->SetBranchAddress("LTr.",&ltr);
2880 Int_t entries= tree->GetEntries();
2881 TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
2882 Double_t bz=AliTrackerBase::GetBz();
2883 //
284418bc 2884 // Double_t globalXYZ[3];
2885 //Double_t globalXYZCorr[3];
46e89793 2886 for (Int_t ientry=0; ientry<entries; ientry++){
2887 tree->GetEntry(ientry);
2888 if (!ltr->GetVecGX()){
2889 ltr->UpdatePoints();
2890 }
2891 //
2892 TVectorD fit10(5);
2893 TVectorD fit5(5);
2894 printf("Entry\t%d\n",ientry);
2895 for (Int_t irow0=0; irow0<158; irow0+=1){
2896 //
2897 TLinearFitter fitter10(4,"hyp3");
2898 TLinearFitter fitter5(2,"hyp1");
2899 Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0];
2900 if (sector<0) continue;
2901 //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue;
2902
2903 Double_t refX= (*ltr->GetVecLX())[irow0];
2904 Int_t firstRow1 = TMath::Max(irow0-10,0);
2905 Int_t lastRow1 = TMath::Min(irow0+10,158);
2906 Double_t padWidth=(irow0<64)?0.4:0.6;
2907 // make long range fit
2908 for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){
2909 if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
2910 if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
2911 if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
2912 Double_t idealX= (*ltr->GetVecLX())[irow1];
2913 Double_t idealY= (*ltr->GetVecLY())[irow1];
284418bc 2914 // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
46e89793 2915 Double_t gx= (*ltr->GetVecGX())[irow1];
2916 Double_t gy= (*ltr->GetVecGY())[irow1];
2917 Double_t gz= (*ltr->GetVecGZ())[irow1];
2918 Double_t measY=(*vecdY)[irow1]+idealY;
2919 Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
2920 // deltaR = R distorted -R ideal
2921 Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
2922 fitter10.AddPoint(xxx,measY,1);
2923 }
2924 Bool_t isOK=kTRUE;
2925 Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
2926 Double_t mean10 =0;// fitter10.GetParameter(0);
2927 Double_t slope10 =0;// fitter10.GetParameter(0);
2928 Double_t cosPart10 = 0;// fitter10.GetParameter(2);
2929 Double_t sinPart10 =0;// fitter10.GetParameter(3);
2930
2931 if (fitter10.GetNpoints()>10){
2932 fitter10.Eval();
2933 rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
2934 mean10 = fitter10.GetParameter(0);
2935 slope10 = fitter10.GetParameter(1);
2936 cosPart10 = fitter10.GetParameter(2);
2937 sinPart10 = fitter10.GetParameter(3);
2938 //
2939 // make short range fit
2940 //
2941 for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){
2942 if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
2943 if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
2944 if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
2945 Double_t idealX= (*ltr->GetVecLX())[irow1];
2946 Double_t idealY= (*ltr->GetVecLY())[irow1];
284418bc 2947 // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
46e89793 2948 Double_t gx= (*ltr->GetVecGX())[irow1];
2949 Double_t gy= (*ltr->GetVecGY())[irow1];
2950 Double_t gz= (*ltr->GetVecGZ())[irow1];
2951 Double_t measY=(*vecdY)[irow1]+idealY;
2952 Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
2953 // deltaR = R distorted -R ideal
2954 Double_t expY= mean10+slope10*(idealX+deltaR-refX);
2955 if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue;
2956 //
2957 Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
2958 Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
2959 fitter5.AddPoint(xxx,measY-corr,1);
2960 }
2961 }else{
2962 isOK=kFALSE;
2963 }
2964 if (fitter5.GetNpoints()<8) isOK=kFALSE;
2965
2966 Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
2967 Double_t offset5 =0;// fitter5.GetParameter(0);
2968 Double_t slope5 =0;// fitter5.GetParameter(0);
2969 if (isOK){
2970 fitter5.Eval();
2971 rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
2972 offset5 = fitter5.GetParameter(0);
2973 slope5 = fitter5.GetParameter(0);
2974 }
2975 //
2976 Double_t dtype=5;
2977 Double_t ptype=0;
2978 Double_t phi =(*ltr->GetVecPhi())[irow0];
2979 Double_t theta =ltr->GetTgl();
2980 Double_t mean=(vecdY)->GetMatrixArray()[irow0];
2981 Double_t gx=0,gy=0,gz=0;
2982 Double_t snp = (*ltr->GetVecP2())[irow0];
2983 Int_t bundle= ltr->GetBundle();
2984 Int_t id= ltr->GetId();
2985 // Double_t rms = err->GetMatrixArray()[irow];
2986 //
2987 gx = (*ltr->GetVecGX())[irow0];
2988 gy = (*ltr->GetVecGY())[irow0];
2989 gz = (*ltr->GetVecGZ())[irow0];
2990 Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0);
2991 fitter10.GetParameters(fit10);
2992 fitter5.GetParameters(fit5);
2993 Double_t idealY= (*ltr->GetVecLY())[irow0];
2994 Double_t measY=(*vecdY)[irow0]+idealY;
2995 Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
2996 if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE;
2997 //
2998 (*pcstream)<<"fitFull"<< // dumpe also intermediate results
2999 "bz="<<bz<< // magnetic filed used
3000 "dtype="<<dtype<< // detector match type
3001 "ptype="<<ptype<< // parameter type
3002 "theta="<<theta<< // theta
3003 "phi="<<phi<< // phi
3004 "snp="<<snp<< // snp
3005 "sector="<<sector<<
3006 "bundle="<<bundle<<
3007// // "dsec="<<dsec<<
3008 "refX="<<refX<< // reference radius
3009 "gx="<<gx<< // global position
3010 "gy="<<gy<< // global position
3011 "gz="<<gz<< // global position
3012 "dRrec="<<dRrec<< // delta Radius in reconstruction
3013 "id="<<id<< //bundle
3014 "rms10="<<rms10<<
3015 "rms5="<<rms5<<
3016 "fit10.="<<&fit10<<
3017 "fit5.="<<&fit5<<
3018 "measY="<<measY<<
3019 "mean="<<mean<<
3020 "idealY="<<idealY<<
3021 "corr="<<corr<<
3022 "isOK="<<isOK<<
3023 "\n";
3024 }
3025 }
3026 delete pcstream;
3027}