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