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