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