1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // _________________________________________________________________
19 // <h2> AliTPCCorrection class </h2>
21 // The AliTPCCorrection class provides a general framework to deal with space point distortions.
22 // An correction class which inherits from here is for example AliTPCExBBShape or AliTPCExBTwist. <br>
23 // General virtual functions are (for example) CorrectPoint(x,roc) where x is the vector of initial
24 // positions in cartesian coordinates and roc represents the read-out chamber number according to
25 // the offline numbering convention. The vector x is overwritten with the corrected coordinates. <br>
26 // An alternative usage would be CorrectPoint(x,roc,dx), which leaves the vector x untouched, but
27 // returns the distortions via the vector dx. <br>
28 // This class is normally used via the general class AliTPCComposedCorrection.
30 // Furthermore, the class contains basic geometrical descriptions like field cage radii
31 // (fgkIFCRadius, fgkOFCRadius) and length (fgkTPCZ0) plus the voltages. Also, the definitions
32 // of size and widths of the fulcrums building the grid of the final look-up table, which is
33 // then interpolated, is defined in kNX and fgkXList).
35 // All physics-model classes below are derived from this class in order to not duplicate code
36 // and to allow a uniform treatment of all physics models.
38 // <h3> Poisson solver </h3>
39 // A numerical solver of the Poisson equation (relaxation technique) is implemented for 2-dimensional
40 // geometries (r,z) as well as for 3-dimensional problems (r,$\phi$,z). The corresponding function
41 // names are PoissonRelaxation?D. The relevant function arguments are the arrays of the boundary and
42 // initial conditions (ArrayofArrayV, ArrayofChargeDensities) as well as the grid granularity which
43 // is used during the calculation. These inputs can be chosen according to the needs of the physical
44 // effect which is supposed to be simulated. In the 3D version, different symmetry conditions can be set
45 // in order to reduce the calculation time (used in AliTPCFCVoltError3D).
47 // <h3> Unified plotting functionality </h3>
48 // Generic plot functions were implemented. They return a histogram pointer in the chosen plane of
49 // the TPC drift volume with a selectable grid granularity and the magnitude of the correction vector.
50 // For example, the function CreateHistoDZinXY(z,nx,ny) returns a 2-dimensional histogram which contains
51 // the longitudinal corrections $dz$ in the (x,y)-plane at the given z position with the granularity of
52 // nx and ny. The magnitude of the corrections is defined by the class from which this function is called.
53 // In the same manner, standard plots for the (r,$\phi$)-plane and for the other corrections like $dr$ and $rd\phi$ are available
55 // Note: This class is normally used via the class AliTPCComposedCorrection
58 // Begin_Macro(source)
60 // gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
61 // TCanvas *c2 = new TCanvas("cAliTPCCorrection","cAliTPCCorrection",700,1050); c2->Divide(2,3);
62 // AliTPCROCVoltError3D roc; // EXAMPLE PLOTS - SEE BELOW
63 // roc.SetROCDataFileName("$ALICE_ROOT/TPC/Calib/maps/TPCROCdzSurvey.root");
64 // roc.SetOmegaTauT1T2(0,1,1); // B=0
65 // Float_t z0 = 1; // at +1 cm -> A side
66 // c2->cd(1); roc.CreateHistoDRinXY(1.,300,300)->Draw("cont4z");
67 // c2->cd(3);roc.CreateHistoDRPhiinXY(1.,300,300)->Draw("cont4z");
68 // c2->cd(5);roc.CreateHistoDZinXY(1.,300,300)->Draw("cont4z");
70 // c2->cd(2);roc.CreateHistoDRinZR(phi0)->Draw("surf2");
71 // c2->cd(4);roc.CreateHistoDRPhiinZR(phi0)->Draw("surf2");
72 // c2->cd(6);roc.CreateHistoDZinZR(phi0)->Draw("surf2");
79 // Date: 27/04/2010 <br>
80 // Authors: Magnus Mager, Stefan Rossegger, Jim Thomas
82 // _________________________________________________________________
85 #include "Riostream.h"
90 #include <TTreeStream.h>
93 #include <TTimeStamp.h>
94 #include <AliCDBStorage.h>
96 #include <AliCDBMetaData.h>
98 #include "AliTPCParamSR.h"
100 #include "AliTPCCorrection.h"
103 #include "AliExternalTrackParam.h"
104 #include "AliTrackPointArray.h"
105 #include "TDatabasePDG.h"
106 #include "AliTrackerBase.h"
107 #include "AliTPCROC.h"
108 #include "THnSparse.h"
110 #include "AliTPCLaserTrack.h"
111 #include "AliESDVertex.h"
112 #include "AliVertexerTracks.h"
113 #include "TDatabasePDG.h"
117 #include "TDatabasePDG.h"
119 #include "AliTPCTransform.h"
120 #include "AliTPCcalibDB.h"
121 #include "AliTPCExB.h"
123 //#include "AliTPCRecoParam.h"
124 #include "TLinearFitter.h"
125 #include <AliSysInfo.h>
127 ClassImp(AliTPCCorrection)
130 TObjArray *AliTPCCorrection::fgVisualCorrection=0;
131 // instance of correction for visualization
134 // FIXME: the following values should come from the database
135 const Double_t AliTPCCorrection::fgkTPCZ0 = 249.7; // nominal gating grid position
136 const Double_t AliTPCCorrection::fgkIFCRadius= 83.5; // radius which renders the "18 rod manifold" best -> compare calc. of Jim Thomas
137 // compare gkIFCRadius= 83.05: Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
138 const Double_t AliTPCCorrection::fgkOFCRadius= 254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
139 const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
140 const Double_t AliTPCCorrection::fgkCathodeV = -100000.0; // Cathode Voltage (volts)
141 const Double_t AliTPCCorrection::fgkGG = -70.0; // Gating Grid voltage (volts)
143 const Double_t AliTPCCorrection::fgkdvdE = 0.0024; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
145 const Double_t AliTPCCorrection::fgkEM = -1.602176487e-19/9.10938215e-31; // charge/mass in [C/kg]
146 const Double_t AliTPCCorrection::fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)]
149 AliTPCCorrection::AliTPCCorrection()
150 : TNamed("correction_unity","unity"),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1), fIsLocal(kFALSE)
153 // default constructor
155 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
157 InitLookUpfulcrums();
161 AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
162 : TNamed(name,title),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1), fIsLocal(kFALSE)
165 // default constructor, that set the name and title
167 if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
169 InitLookUpfulcrums();
173 AliTPCCorrection::~AliTPCCorrection() {
175 // virtual destructor
179 Bool_t AliTPCCorrection::AddCorrectionCompact(AliTPCCorrection* corr, Double_t weight){
181 // Add correction and make them compact
183 // - origin of distortion/correction are additive
184 // - only correction ot the same type supported ()
186 AliError("Zerro pointer - correction");
189 AliError(TString::Format("Correction %s not implementend",IsA()->GetName()).Data());
194 void AliTPCCorrection::CorrectPoint(Float_t x[], Short_t roc) {
196 // Corrects the initial coordinates x (cartesian coordinates)
197 // according to the given effect (inherited classes)
198 // roc represents the TPC read out chamber (offline numbering convention)
201 GetCorrection(x,roc,dx);
202 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
205 void AliTPCCorrection::CorrectPoint(const Float_t x[], Short_t roc,Float_t xp[]) {
207 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
208 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
209 // roc represents the TPC read out chamber (offline numbering convention)
212 GetCorrection(x,roc,dx);
213 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
216 void AliTPCCorrection::DistortPoint(Float_t x[], Short_t roc) {
218 // Distorts the initial coordinates x (cartesian coordinates)
219 // according to the given effect (inherited classes)
220 // roc represents the TPC read out chamber (offline numbering convention)
223 GetDistortion(x,roc,dx);
224 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
227 void AliTPCCorrection::DistortPointLocal(Float_t x[], Short_t roc) {
229 // Distorts the initial coordinates x (cartesian coordinates)
230 // according to the given effect (inherited classes)
231 // roc represents the TPC read out chamber (offline numbering convention)
233 Float_t gxyz[3]={0,0,0};
234 Double_t alpha = TMath::TwoPi()*(roc%18+0.5)/18;
235 Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
236 gxyz[0]= ca*x[0]+sa*x[1];
237 gxyz[1]= -sa*x[0]+ca*x[1];
239 DistortPoint(gxyz,roc);
240 x[0]= ca*gxyz[0]-sa*gxyz[1];
241 x[1]= +sa*gxyz[0]+ca*gxyz[1];
244 void AliTPCCorrection::CorrectPointLocal(Float_t x[], Short_t roc) {
246 // Distorts the initial coordinates x (cartesian coordinates)
247 // according to the given effect (inherited classes)
248 // roc represents the TPC read out chamber (offline numbering convention)
250 Float_t gxyz[3]={0,0,0};
251 Double_t alpha = TMath::TwoPi()*(roc%18+0.5)/18;
252 Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
253 gxyz[0]= ca*x[0]+sa*x[1];
254 gxyz[1]= -sa*x[0]+ca*x[1];
256 CorrectPoint(gxyz,roc);
257 x[0]= ca*gxyz[0]-sa*gxyz[1];
258 x[1]= sa*gxyz[0]+ca*gxyz[1];
262 void AliTPCCorrection::DistortPoint(const Float_t x[], Short_t roc,Float_t xp[]) {
264 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
265 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
266 // roc represents the TPC read out chamber (offline numbering convention)
269 GetDistortion(x,roc,dx);
270 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
273 void AliTPCCorrection::GetCorrection(const Float_t /*x*/[], Short_t /*roc*/,Float_t dx[]) {
275 // This function delivers the correction values dx in respect to the inital coordinates x
276 // roc represents the TPC read out chamber (offline numbering convention)
277 // Note: The dx is overwritten by the inherited effectice class ...
279 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
282 void AliTPCCorrection::GetDistortion(const Float_t x[], Short_t roc,Float_t dx[]) {
284 // This function delivers the distortion values dx in respect to the inital coordinates x
285 // roc represents the TPC read out chamber (offline numbering convention)
287 GetCorrection(x,roc,dx);
288 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
291 void AliTPCCorrection::GetCorrectionDz(const Float_t x[], Short_t roc,Float_t dx[], Float_t delta) {
292 // author: marian.ivanov@cern.ch
294 // In this (virtual)function calculates the dx'/dz, dy'/dz and dz'/dz at given point (x,y,z)
295 // Generic implementation. Better precision can be acchieved knowing the internal structure
296 // of underlying trasnformation. Derived classes can reimplement it.
297 // To calculate correction is fitted in small neighberhood:
298 // (x+-delta,y+-delta,z+-delta) where delta is an argument
301 // x[] - space point corrdinate
302 // roc - readout chamber identifier (important e.g to do not miss the side of detector)
303 // delta - define the size of neighberhood
305 // dx[] - array {dx'/dz, dy'/dz , dz'/dz }
307 // if (fIsLocal){ //standard implemenation provides the correction/distortion integrated over full drift length
310 // GetCorrection(xyz,roc,dxyz);
312 static TLinearFitter fitx(2,"pol1");
313 static TLinearFitter fity(2,"pol1");
314 static TLinearFitter fitz(2,"pol1");
320 //adjust limits around CE to stay on one side
323 if ((x[2]+zmin*delta)<0){
335 if ((x[2]+zmax*delta)>0){
345 for (Int_t xdelta=-1; xdelta<=1; xdelta++)
346 for (Int_t ydelta=-1; ydelta<=1; ydelta++){
347 // for (Int_t zdelta=-1; zdelta<=1; zdelta++){
348 // for (Int_t xdelta=-2; xdelta<=0; xdelta++)
349 // for (Int_t ydelta=-2; ydelta<=0; ydelta++){
350 for (Int_t zdelta=zmin; zdelta<=zmax; zdelta++){
351 //TODO: what happens if x[2] is on the A-Side, but x[2]+zdelta*delta
352 // will be on the C-Side?
353 Float_t xyz[3]={x[0]+xdelta*delta, x[1]+ydelta*delta, x[2]+zdelta*delta};
355 GetCorrection(xyz,roc,dxyz);
356 Double_t adelta=zdelta*delta;
357 fitx.AddPoint(&adelta, dxyz[0]);
358 fity.AddPoint(&adelta, dxyz[1]);
359 fitz.AddPoint(&adelta, dxyz[2]);
365 dx[0] = fitx.GetParameter(1);
366 dx[1] = fity.GetParameter(1);
367 dx[2] = fitz.GetParameter(1);
370 void AliTPCCorrection::GetDistortionDz(const Float_t x[], Short_t roc,Float_t dx[], Float_t delta) {
371 // author: marian.ivanov@cern.ch
373 // In this (virtual)function calculates the dx'/dz, dy'/dz and dz'/dz at given point (x,y,z)
374 // Generic implementation. Better precision can be acchieved knowing the internal structure
375 // of underlying trasnformation. Derived classes can reimplement it.
376 // To calculate distortion is fitted in small neighberhood:
377 // (x+-delta,y+-delta,z+-delta) where delta is an argument
380 // x[] - space point corrdinate
381 // roc - readout chamber identifier (important e.g to do not miss the side of detector)
382 // delta - define the size of neighberhood
384 // dx[] - array {dx'/dz, dy'/dz , dz'/dz }
386 static TLinearFitter fitx(2,"pol1");
387 static TLinearFitter fity(2,"pol1");
388 static TLinearFitter fitz(2,"pol1");
395 //adjust limits around CE to stay on one side
398 if ((x[2]+zmin*delta)<0){
404 if ((x[2]+zmax*delta)>0){
410 //TODO: in principle one shuld check that x[2]+zdelta*delta does not get 'out of' bounds,
411 // so close to the CE it doesn't change the sign, since then the corrections will be wrong ...
412 for (Int_t xdelta=-1; xdelta<=1; xdelta++)
413 for (Int_t ydelta=-1; ydelta<=1; ydelta++){
414 for (Int_t zdelta=zmin; zdelta<=zmax; zdelta++){
415 //TODO: what happens if x[2] is on the A-Side, but x[2]+zdelta*delta
416 // will be on the C-Side?
417 //TODO: For the C-Side, does this have the correct sign?
418 Float_t xyz[3]={x[0]+xdelta*delta, x[1]+ydelta*delta, x[2]+zdelta*delta};
420 GetDistortion(xyz,roc,dxyz);
421 Double_t adelta=zdelta*delta;
422 fitx.AddPoint(&adelta, dxyz[0]);
423 fity.AddPoint(&adelta, dxyz[1]);
424 fitz.AddPoint(&adelta, dxyz[2]);
430 dx[0] = fitx.GetParameter(1);
431 dx[1] = fity.GetParameter(1);
432 dx[2] = fitz.GetParameter(1);
435 void AliTPCCorrection::GetCorrectionIntegralDz(const Float_t x[], Short_t roc,Float_t dx[], Float_t delta){
437 // Integrate 3D distortion along drift lines starting from the roc plane
438 // to the expected z position of the point, this assumes that dz is small
439 // and the error propagating to z' instead of the correct z is negligible
440 // To define the drift lines virtual function AliTPCCorrection::GetCorrectionDz is used
443 // x[] - space point corrdinate
444 // roc - readout chamber identifier (important e.g to do not miss the side of detector)
445 // delta - define the size of neighberhood
447 // dx[] - array { integral(dx'/dz), integral(dy'/dz) , integral(dz'/dz) }
449 Float_t zroc= ((roc%36)<18) ? fgkTPCZ0:-fgkTPCZ0;
450 Double_t zdrift = TMath::Abs(x[2]-zroc);
451 Int_t nsteps = Int_t(zdrift/delta)+1;
454 Float_t xyz[3]={x[0],x[1],zroc};
455 Float_t dxyz[3]={x[0],x[1],x[2]};
456 Short_t side=(roc/18)%2;
457 Float_t sign=1-2*side;
459 for (Int_t i=0;i<nsteps; i++){
460 //propagate backwards, therefore opposite signs
461 Float_t deltaZ=delta*(-sign);
462 // if (xyz[2]+deltaZ>fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-fgkTPCZ0);
463 // if (xyz[2]-deltaZ<-fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-fgkTPCZ0);
464 // protect again integrating through the CE
466 if (xyz[2]+deltaZ<0) deltaZ=-xyz[2]+1e-20;
468 if (xyz[2]+deltaZ>0) deltaZ=xyz[2]-+1e-20;
470 // since at larger drift (smaller z) the corrections are larger (absolute, but negative)
471 // the slopes will be positive.
472 // but since we chose deltaZ opposite sign the singn of the corretion should be fine
474 Float_t xyz2[3]={xyz[0],xyz[1],static_cast<Float_t>(xyz[2]+deltaZ/2.)};
475 GetCorrectionDz(xyz2,roc,dxyz,delta/2.);
476 xyz[0]+=deltaZ*dxyz[0];
477 xyz[1]+=deltaZ*dxyz[1];
479 sumdz+=deltaZ*dxyz[2];
484 dx[2]= sumdz; //TODO: is sumdz correct?
487 void AliTPCCorrection::GetDistortionIntegralDz(const Float_t x[], Short_t roc,Float_t dx[], Float_t delta){
489 // Integrate 3D distortion along drift lines
490 // To define the drift lines virtual function AliTPCCorrection::GetCorrectionDz is used
493 // x[] - space point corrdinate
494 // roc - readout chamber identifier (important e.g to do not miss the side of detector)
495 // delta - define the size of neighberhood
497 // dx[] - array { integral(dx'/dz), integral(dy'/dz) , integral(dz'/dz) }
499 Float_t zroc= ((roc%36)<18) ? fgkTPCZ0:-fgkTPCZ0;
500 Double_t zdrift = TMath::Abs(x[2]-zroc);
501 Int_t nsteps = Int_t(zdrift/delta)+1;
504 Float_t xyz[3]={x[0],x[1],x[2]};
505 Float_t dxyz[3]={x[0],x[1],x[2]};
506 Float_t sign=((roc%36)<18) ? 1.:-1.;
508 for (Int_t i=0;i<nsteps; i++){
509 Float_t deltaZ=delta;
510 if (xyz[2]+deltaZ>fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-zroc);
511 if (xyz[2]-deltaZ<-fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-zroc);
512 // since at larger drift (smaller z) the distortions are larger
513 // the slopes will be negative.
514 // and since we are moving towards the read-out plane the deltaZ for
515 // weighting the dK/dz should have the opposite sign
517 Float_t xyz2[3]={xyz[0],xyz[1],static_cast<Float_t>(xyz[2]+deltaZ/2.)};
518 GetDistortionDz(xyz2,roc,dxyz,delta/2.);
519 xyz[0]+=-deltaZ*dxyz[0];
520 xyz[1]+=-deltaZ*dxyz[1];
521 xyz[2]+=deltaZ; //TODO: Should this also be corrected for the dxyz[2]
522 sumdz+=-deltaZ*dxyz[2];
527 dx[2]= sumdz; //TODO: is sumdz correct?
532 void AliTPCCorrection::Init() {
534 // Initialization funtion (not used at the moment)
538 void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
544 void AliTPCCorrection::Print(Option_t* /*option*/) const {
546 // Print function to check which correction classes are used
547 // option=="d" prints details regarding the setted magnitude
548 // option=="a" prints the C0 and C1 coefficents for calibration purposes
550 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
553 void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
555 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
556 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
561 //SetOmegaTauT1T2(omegaTau, t1, t2);
564 TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
566 // Simple plot functionality.
567 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
568 // in respect to position z within the XY plane.
569 // The histogramm has nx times ny entries.
571 AliTPCParam* tpcparam = new AliTPCParamSR;
573 TH2F *h=CreateTH2F("dr_xy", TString::Format("%s: DRinXY Z=%2.0f", GetTitle(),z).Data(),"x [cm]","y [cm]","dr [cm]",
574 nx,-250.,250.,ny,-250.,250.);
577 Int_t roc=z>0.?0:18; // FIXME
578 for (Int_t iy=1;iy<=ny;++iy) {
579 x[1]=h->GetYaxis()->GetBinCenter(iy);
580 for (Int_t ix=1;ix<=nx;++ix) {
581 x[0]=h->GetXaxis()->GetBinCenter(ix);
582 GetCorrection(x,roc,dx);
583 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
584 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
585 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
586 h->SetBinContent(ix,iy,r1-r0);
589 h->SetBinContent(ix,iy,0.);
596 TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
598 // Simple plot functionality.
599 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
600 // in respect to position z within the XY plane.
601 // The histogramm has nx times ny entries.
604 AliTPCParam* tpcparam = new AliTPCParamSR;
606 TH2F *h=CreateTH2F("drphi_xy",TString::Format("%s: DRPhiinXY Z=%2.0f", GetTitle(),z).Data(),"x [cm]","y [cm]","drphi [cm]",
607 nx,-250.,250.,ny,-250.,250.);
610 Int_t roc=z>0.?0:18; // FIXME
611 for (Int_t iy=1;iy<=ny;++iy) {
612 x[1]=h->GetYaxis()->GetBinCenter(iy);
613 for (Int_t ix=1;ix<=nx;++ix) {
614 x[0]=h->GetXaxis()->GetBinCenter(ix);
615 GetCorrection(x,roc,dx);
616 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
617 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
618 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
619 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
621 Float_t dphi=phi1-phi0;
622 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
623 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
625 h->SetBinContent(ix,iy,r0*dphi);
628 h->SetBinContent(ix,iy,0.);
635 TH2F* AliTPCCorrection::CreateHistoDZinXY(Float_t z,Int_t nx,Int_t ny) {
637 // Simple plot functionality.
638 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
639 // in respect to position z within the XY plane.
640 // The histogramm has nx times ny entries.
643 AliTPCParam* tpcparam = new AliTPCParamSR;
645 TH2F *h=CreateTH2F("dz_xy",TString::Format("%s: DZinXY Z=%2.0f", GetTitle(),z).Data(),"x [cm]","y [cm]","dz [cm]",
646 nx,-250.,250.,ny,-250.,250.);
649 Int_t roc=z>0.?0:18; // FIXME
650 for (Int_t iy=1;iy<=ny;++iy) {
651 x[1]=h->GetYaxis()->GetBinCenter(iy);
652 for (Int_t ix=1;ix<=nx;++ix) {
653 x[0]=h->GetXaxis()->GetBinCenter(ix);
654 GetCorrection(x,roc,dx);
655 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
656 if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
657 h->SetBinContent(ix,iy,dx[2]);
660 h->SetBinContent(ix,iy,0.);
667 TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
669 // Simple plot functionality.
670 // Returns a 2d hisogram which represents the corrections in r direction (dr)
671 // in respect to angle phi within the ZR plane.
672 // The histogramm has nx times ny entries.
674 TH2F *h=CreateTH2F("dr_zr",TString::Format("%s: DRinZR Phi=%2.2f", GetTitle(),phi).Data(),"z [cm]","r [cm]","dr [cm]",
675 nz,-250.,250.,nr,85.,250.);
677 for (Int_t ir=1;ir<=nr;++ir) {
678 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
679 x[0]=radius*TMath::Cos(phi);
680 x[1]=radius*TMath::Sin(phi);
681 for (Int_t iz=1;iz<=nz;++iz) {
682 x[2]=h->GetXaxis()->GetBinCenter(iz);
683 Int_t roc=x[2]>0.?0:18; // FIXME
684 GetCorrection(x,roc,dx);
685 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
686 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
687 h->SetBinContent(iz,ir,r1-r0);
694 TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
696 // Simple plot functionality.
697 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
698 // in respect to angle phi within the ZR plane.
699 // The histogramm has nx times ny entries.
701 TH2F *h=CreateTH2F("drphi_zr", TString::Format("%s: DRPhiinZR R=%2.2f", GetTitle(),phi).Data(),"z [cm]","r [cm]","drphi [cm]",
702 nz,-250.,250.,nr,85.,250.);
704 for (Int_t iz=1;iz<=nz;++iz) {
705 x[2]=h->GetXaxis()->GetBinCenter(iz);
706 Int_t roc=x[2]>0.?0:18; // FIXME
707 for (Int_t ir=1;ir<=nr;++ir) {
708 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
709 x[0]=radius*TMath::Cos(phi);
710 x[1]=radius*TMath::Sin(phi);
711 GetCorrection(x,roc,dx);
712 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
713 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
714 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
716 Float_t dphi=phi1-phi0;
717 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
718 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
720 h->SetBinContent(iz,ir,r0*dphi);
726 TH2F* AliTPCCorrection::CreateHistoDZinZR(Float_t phi,Int_t nz,Int_t nr) {
728 // Simple plot functionality.
729 // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
730 // in respect to angle phi within the ZR plane.
731 // The histogramm has nx times ny entries.
733 TH2F *h=CreateTH2F("dz_zr",TString::Format("%s: DZinZR Z=%2.0f", GetTitle(),phi).Data(),"z [cm]","r [cm]","dz [cm]",
734 nz,-250.,250.,nr,85.,250.);
736 for (Int_t ir=1;ir<=nr;++ir) {
737 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
738 x[0]=radius*TMath::Cos(phi);
739 x[1]=radius*TMath::Sin(phi);
740 for (Int_t iz=1;iz<=nz;++iz) {
741 x[2]=h->GetXaxis()->GetBinCenter(iz);
742 Int_t roc=x[2]>0.?0:18; // FIXME
743 GetCorrection(x,roc,dx);
744 h->SetBinContent(iz,ir,dx[2]);
752 TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
753 const char *xlabel,const char *ylabel,const char *zlabel,
754 Int_t nbinsx,Double_t xlow,Double_t xup,
755 Int_t nbinsy,Double_t ylow,Double_t yup) {
757 // Helper function to create a 2d histogramm of given size
763 while (gDirectory->FindObject(hname.Data())) {
770 TH2F *h=new TH2F(hname.Data(),title,
773 h->GetXaxis()->SetTitle(xlabel);
774 h->GetYaxis()->SetTitle(ylabel);
775 h->GetZaxis()->SetTitle(zlabel);
780 // Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
782 void AliTPCCorrection::Interpolate2DEdistortion( Int_t order, Double_t r, Double_t z,
783 const Double_t er[kNZ][kNR], Double_t &erValue ) {
785 // Interpolate table - 2D interpolation
787 Double_t saveEr[5] = {0,0,0,0,0};
789 Search( kNZ, fgkZList, z, fJLow ) ;
790 Search( kNR, fgkRList, r, fKLow ) ;
791 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
792 if ( fKLow < 0 ) fKLow = 0 ;
793 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
794 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
796 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
797 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
799 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
803 void AliTPCCorrection::Interpolate3DEdistortion( Int_t order, Double_t r, Float_t phi, Double_t z,
804 const Double_t er[kNZ][kNPhi][kNR], const Double_t ephi[kNZ][kNPhi][kNR], const Double_t ez[kNZ][kNPhi][kNR],
805 Double_t &erValue, Double_t &ephiValue, Double_t &ezValue) {
807 // Interpolate table - 3D interpolation
810 Double_t saveEr[5]= {0,0,0,0,0};
811 Double_t savedEr[5]= {0,0,0,0,0} ;
813 Double_t saveEphi[5]= {0,0,0,0,0};
814 Double_t savedEphi[5]= {0,0,0,0,0} ;
816 Double_t saveEz[5]= {0,0,0,0,0};
817 Double_t savedEz[5]= {0,0,0,0,0} ;
819 Search( kNZ, fgkZList, z, fILow ) ;
820 Search( kNPhi, fgkPhiList, z, fJLow ) ;
821 Search( kNR, fgkRList, r, fKLow ) ;
823 if ( fILow < 0 ) fILow = 0 ; // check if out of range
824 if ( fJLow < 0 ) fJLow = 0 ;
825 if ( fKLow < 0 ) fKLow = 0 ;
827 if ( fILow + order >= kNZ - 1 ) fILow = kNZ - 1 - order ;
828 if ( fJLow + order >= kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
829 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
831 for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
832 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
833 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r ) ;
834 saveEphi[j-fJLow] = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
835 saveEz[j-fJLow] = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r ) ;
837 savedEr[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi ) ;
838 savedEphi[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
839 savedEz[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi ) ;
841 erValue = Interpolate( &fgkZList[fILow], savedEr, order, z ) ;
842 ephiValue = Interpolate( &fgkZList[fILow], savedEphi, order, z ) ;
843 ezValue = Interpolate( &fgkZList[fILow], savedEz, order, z ) ;
847 Double_t AliTPCCorrection::Interpolate2DTable( Int_t order, Double_t x, Double_t y,
848 Int_t nx, Int_t ny, const Double_t xv[], const Double_t yv[],
849 const TMatrixD &array ) {
851 // Interpolate table (TMatrix format) - 2D interpolation
854 static Int_t jlow = 0, klow = 0 ;
855 Double_t saveArray[5] = {0,0,0,0,0} ;
857 Search( nx, xv, x, jlow ) ;
858 Search( ny, yv, y, klow ) ;
859 if ( jlow < 0 ) jlow = 0 ; // check if out of range
860 if ( klow < 0 ) klow = 0 ;
861 if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
862 if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
864 for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
866 Double_t *ajkl = &((TMatrixD&)array)(j,klow);
867 saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
870 return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
874 Double_t AliTPCCorrection::Interpolate3DTable( Int_t order, Double_t x, Double_t y, Double_t z,
875 Int_t nx, Int_t ny, Int_t nz,
876 const Double_t xv[], const Double_t yv[], const Double_t zv[],
877 TMatrixD **arrayofArrays ) {
879 // Interpolate table (TMatrix format) - 3D interpolation
882 static Int_t ilow = 0, jlow = 0, klow = 0 ;
883 Double_t saveArray[5]= {0,0,0,0,0};
884 Double_t savedArray[5]= {0,0,0,0,0} ;
886 Search( nx, xv, x, ilow ) ;
887 Search( ny, yv, y, jlow ) ;
888 Search( nz, zv, z, klow ) ;
890 if ( ilow < 0 ) ilow = 0 ; // check if out of range
891 if ( jlow < 0 ) jlow = 0 ;
892 if ( klow < 0 ) klow = 0 ;
894 if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
895 if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
896 if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
898 for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
900 TMatrixD &table = *arrayofArrays[k] ;
901 for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
903 saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
905 savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
907 return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
911 Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
912 Int_t order, Double_t x ) {
914 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
918 if ( order == 2 ) { // Quadratic Interpolation = 2
919 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
920 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
921 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
922 } else { // Linear Interpolation = 1
923 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
930 Float_t AliTPCCorrection::Interpolate2DTable( Int_t order, Double_t x, Double_t y,
931 Int_t nx, Int_t ny, const Double_t xv[], const Double_t yv[],
932 const TMatrixF &array ) {
934 // Interpolate table (TMatrix format) - 2D interpolation
935 // Float version (in order to decrease the OCDB size)
938 static Int_t jlow = 0, klow = 0 ;
939 Float_t saveArray[5] = {0.,0.,0.,0.,0.} ;
941 Search( nx, xv, x, jlow ) ;
942 Search( ny, yv, y, klow ) ;
943 if ( jlow < 0 ) jlow = 0 ; // check if out of range
944 if ( klow < 0 ) klow = 0 ;
945 if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
946 if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
948 for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
950 Float_t *ajkl = &((TMatrixF&)array)(j,klow);
951 saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
954 return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
958 Float_t AliTPCCorrection::Interpolate3DTable( Int_t order, Double_t x, Double_t y, Double_t z,
959 Int_t nx, Int_t ny, Int_t nz,
960 const Double_t xv[], const Double_t yv[], const Double_t zv[],
961 TMatrixF **arrayofArrays ) {
963 // Interpolate table (TMatrix format) - 3D interpolation
964 // Float version (in order to decrease the OCDB size)
967 static Int_t ilow = 0, jlow = 0, klow = 0 ;
968 Float_t saveArray[5]= {0.,0.,0.,0.,0.};
969 Float_t savedArray[5]= {0.,0.,0.,0.,0.} ;
971 Search( nx, xv, x, ilow ) ;
972 Search( ny, yv, y, jlow ) ;
973 Search( nz, zv, z, klow ) ;
975 if ( ilow < 0 ) ilow = 0 ; // check if out of range
976 if ( jlow < 0 ) jlow = 0 ;
977 if ( klow < 0 ) klow = 0 ;
979 if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
980 if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
981 if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
983 for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
985 TMatrixF &table = *arrayofArrays[k] ;
986 for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
988 saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
990 savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
992 return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
995 Float_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Float_t yArray[],
996 Int_t order, Double_t x ) {
998 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
999 // Float version (in order to decrease the OCDB size)
1003 if ( order == 2 ) { // Quadratic Interpolation = 2
1004 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
1005 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
1006 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
1007 } else { // Linear Interpolation = 1
1008 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
1017 void AliTPCCorrection::Search( Int_t n, const Double_t xArray[], Double_t x, Int_t &low ) {
1019 // Search an ordered table by starting at the most recently used point
1022 Long_t middle, high ;
1023 Int_t ascend = 0, increment = 1 ;
1025 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
1027 if ( low < 0 || low > n-1 ) {
1028 low = -1 ; high = n ;
1029 } else { // Ordered Search phase
1030 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
1031 if ( low == n-1 ) return ;
1033 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
1036 high = low + increment ;
1037 if ( high > n-1 ) { high = n ; break ; }
1040 if ( low == 0 ) { low = -1 ; return ; }
1042 while ( (Int_t)( x < xArray[low] ) == ascend ) {
1045 if ( increment >= high ) { low = -1 ; break ; }
1046 else low = high - increment ;
1051 while ( (high-low) != 1 ) { // Binary Search Phase
1052 middle = ( high + low ) / 2 ;
1053 if ( (Int_t)( x >= xArray[middle] ) == ascend )
1059 if ( x == xArray[n-1] ) low = n-2 ;
1060 if ( x == xArray[0] ) low = 0 ;
1064 void AliTPCCorrection::InitLookUpfulcrums() {
1066 // Initialization of interpolation points - for main look up table
1067 // (course grid in the middle, fine grid on the borders)
1070 AliTPCROC * roc = AliTPCROC::Instance();
1071 const Double_t rLow = TMath::Floor(roc->GetPadRowRadii(0,0))-1; // first padRow plus some margin
1075 for (Int_t i = 1; i<kNR; i++) {
1076 fgkRList[i] = fgkRList[i-1] + 3.5; // 3.5 cm spacing
1077 if (fgkRList[i]<90 ||fgkRList[i]>245)
1078 fgkRList[i] = fgkRList[i-1] + 0.5; // 0.5 cm spacing
1079 else if (fgkRList[i]<100 || fgkRList[i]>235)
1080 fgkRList[i] = fgkRList[i-1] + 1.5; // 1.5 cm spacing
1081 else if (fgkRList[i]<120 || fgkRList[i]>225)
1082 fgkRList[i] = fgkRList[i-1] + 2.5; // 2.5 cm spacing
1086 fgkZList[0] = -249.5;
1087 fgkZList[kNZ-1] = 249.5;
1088 for (Int_t j = 1; j<kNZ/2; j++) {
1089 fgkZList[j] = fgkZList[j-1];
1090 if (TMath::Abs(fgkZList[j])< 0.15)
1091 fgkZList[j] = fgkZList[j-1] + 0.09; // 0.09 cm spacing
1092 else if(TMath::Abs(fgkZList[j])< 0.6)
1093 fgkZList[j] = fgkZList[j-1] + 0.4; // 0.4 cm spacing
1094 else if (TMath::Abs(fgkZList[j])< 2.5 || TMath::Abs(fgkZList[j])>248)
1095 fgkZList[j] = fgkZList[j-1] + 0.5; // 0.5 cm spacing
1096 else if (TMath::Abs(fgkZList[j])<10 || TMath::Abs(fgkZList[j])>235)
1097 fgkZList[j] = fgkZList[j-1] + 1.5; // 1.5 cm spacing
1098 else if (TMath::Abs(fgkZList[j])<25 || TMath::Abs(fgkZList[j])>225)
1099 fgkZList[j] = fgkZList[j-1] + 2.5; // 2.5 cm spacing
1101 fgkZList[j] = fgkZList[j-1] + 4; // 4 cm spacing
1103 fgkZList[kNZ-j-1] = -fgkZList[j];
1107 for (Int_t k = 0; k<kNPhi; k++)
1108 fgkPhiList[k] = TMath::TwoPi()*k/(kNPhi-1);
1114 void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
1115 TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
1116 Int_t rows, Int_t columns, Int_t iterations,
1117 Bool_t rocDisplacement ) {
1119 // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
1121 // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
1122 // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
1123 // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
1124 // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
1125 // you wish to solve Laplaces equation however it should not contain random numbers or you will get
1126 // random numbers back as a solution.
1127 // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
1128 // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
1129 // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
1130 // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
1131 // number of points and solves the problem again. This happens several times until the maximum number
1132 // of points has been included in the array.
1134 // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
1135 // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
1137 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
1139 // Original code by Jim Thomas (STAR TPC Collaboration)
1142 Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
1144 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
1145 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
1146 const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
1148 TMatrixD arrayEr(rows,columns) ;
1149 TMatrixD arrayEz(rows,columns) ;
1151 //Check that number of rows and columns is suitable for a binary expansion
1153 if ( !IsPowerOfTwo(rows-1) ) {
1154 AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
1157 if ( !IsPowerOfTwo(columns-1) ) {
1158 AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
1162 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
1163 // Allow for different size grid spacing in R and Z directions
1164 // Use a binary expansion of the size of the matrix to speed up the solution of the problem
1166 Int_t iOne = (rows-1)/4 ;
1167 Int_t jOne = (columns-1)/4 ;
1168 // Solve for N in 2**N, add one.
1169 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
1171 for ( Int_t count = 0 ; count < loops ; count++ ) {
1172 // Loop while the matrix expands & the resolution increases.
1174 Float_t tempGridSizeR = gridSizeR * iOne ;
1175 Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
1176 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
1178 // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1179 std::vector<float> coef1(rows) ;
1180 std::vector<float> coef2(rows) ;
1182 for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
1183 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1184 coef1[i] = 1.0 + tempGridSizeR/(2*radius);
1185 coef2[i] = 1.0 - tempGridSizeR/(2*radius);
1188 TMatrixD sumChargeDensity(rows,columns) ;
1190 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
1191 Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
1192 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
1193 if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1195 // Add up all enclosed charge density contributions within 1/2 unit in all directions
1196 Float_t weight = 0.0 ;
1198 sumChargeDensity(i,j) = 0.0 ;
1199 for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
1200 for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
1201 if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
1204 // Note that this is cylindrical geometry
1205 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1206 sum += weight*radius ;
1209 sumChargeDensity(i,j) /= sum ;
1211 sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
1215 for ( Int_t k = 1 ; k <= iterations; k++ ) {
1216 // Solve Poisson's Equation
1217 // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
1218 // as interations increase.
1219 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1220 Float_t overRelaxM1 = overRelax - 1.0 ;
1221 Float_t overRelaxtempFourth, overRelaxcoef5 ;
1222 overRelaxtempFourth = overRelax * tempFourth ;
1223 overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
1225 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
1226 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
1228 arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
1229 + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
1230 - overRelaxcoef5 * arrayV(i,j)
1231 + coef1[i] * arrayV(i+iOne,j)
1232 + sumChargeDensity(i,j)
1233 ) * overRelaxtempFourth;
1237 if ( k == iterations ) {
1238 // After full solution is achieved, copy low resolution solution into higher res array
1239 for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
1240 for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
1243 arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
1244 if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
1247 arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
1248 if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
1250 if ( iOne > 1 && jOne > 1 ) {
1251 arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
1252 if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
1253 if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
1254 // Note that this leaves a point at the upper left and lower right corners uninitialized.
1255 // -> Not a big deal.
1264 iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
1265 jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
1267 sumChargeDensity.Clear();
1270 // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
1271 for ( Int_t j = 0 ; j < columns ; j++ ) {
1272 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1273 arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1274 arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1277 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1278 for ( Int_t i = 0 ; i < rows ; i++) {
1279 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1280 arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1281 arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1284 for ( Int_t i = 0 ; i < rows ; i++) {
1285 // Note: go back and compare to old version of this code. See notes below.
1286 // JT Test ... attempt to divide by real Ez not Ez to first order
1287 for ( Int_t j = 0 ; j < columns ; j++ ) {
1288 arrayEz(i,j) += ezField;
1289 // This adds back the overall Z gradient of the field (main E field component)
1291 // Warning: (-=) assumes you are using an error potetial without the overall Field included
1294 // Integrate Er/Ez from Z to zero
1295 for ( Int_t j = 0 ; j < columns ; j++ ) {
1296 for ( Int_t i = 0 ; i < rows ; i++ ) {
1298 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1299 arrayErOverEz(i,j) = 0.0 ;
1300 arrayDeltaEz(i,j) = 0.0 ;
1302 for ( Int_t k = j ; k < columns ; k++ ) {
1303 arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
1304 arrayDeltaEz(i,j) += index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
1305 if ( index != 4 ) index = 4; else index = 2 ;
1308 arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
1309 arrayDeltaEz(i,j) -= (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
1312 arrayErOverEz(i,j) += (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
1313 -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
1314 arrayDeltaEz(i,j) += (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
1315 -2.5*(arrayEz(i,columns-1)-ezField));
1317 if ( j == columns-2 ) {
1318 arrayErOverEz(i,j) = (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
1319 +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
1320 arrayDeltaEz(i,j) = (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
1321 +1.5*(arrayEz(i,columns-1)-ezField) ) ;
1323 if ( j == columns-1 ) {
1324 arrayErOverEz(i,j) = 0.0 ;
1325 arrayDeltaEz(i,j) = 0.0 ;
1330 // calculate z distortion from the integrated Delta Ez residuals
1331 // and include the aquivalence (Volt to cm) of the ROC shift !!
1333 for ( Int_t j = 0 ; j < columns ; j++ ) {
1334 for ( Int_t i = 0 ; i < rows ; i++ ) {
1336 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1337 arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
1339 // ROC Potential in cm aquivalent
1340 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
1341 if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift; // add the ROC misaligment
1351 void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
1352 TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
1353 Int_t rows, Int_t columns, Int_t phislices,
1354 Float_t deltaphi, Int_t iterations, Int_t symmetry,
1355 Bool_t rocDisplacement ) {
1357 // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
1359 // NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
1360 // The number of rows and COLUMNS can be different.
1363 // COLUMNS == 2**N + 1
1364 // PHISLICES == Arbitrary but greater than 3
1366 // DeltaPhi in Radians
1368 // SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
1369 // = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
1371 // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
1373 const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
1375 const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
1376 const Float_t gridSizePhi = deltaphi ;
1377 const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
1378 const Float_t ratioPhi = gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
1379 const Float_t ratioZ = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
1381 TMatrixD arrayE(rows,columns) ;
1383 // Check that the number of rows and columns is suitable for a binary expansion
1384 if ( !IsPowerOfTwo((rows-1)) ) {
1385 AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
1387 if ( !IsPowerOfTwo((columns-1)) ) {
1388 AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
1390 if ( phislices <= 3 ) {
1391 AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
1393 if ( phislices > 1000 ) {
1394 AliError("Poisson3D phislices > 1000 is not allowed (nor wise) ");
1397 // Solve Poisson's equation in cylindrical coordinates by relaxation technique
1398 // Allow for different size grid spacing in R and Z directions
1399 // Use a binary expansion of the matrix to speed up the solution of the problem
1401 Int_t loops, mplus, mminus, signplus, signminus ;
1402 Int_t ione = (rows-1)/4 ;
1403 Int_t jone = (columns-1)/4 ;
1404 loops = TMath::Max(ione, jone) ; // Calculate the number of loops for the binary expansion
1405 loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ; // Solve for N in 2**N
1407 TMatrixD* arrayofSumChargeDensities[1000] ; // Create temporary arrays to store low resolution charge arrays
1409 for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
1410 AliSysInfo::AddStamp("3DInit", 10,0,0);
1412 for ( Int_t count = 0 ; count < loops ; count++ ) { // START the master loop and do the binary expansion
1413 AliSysInfo::AddStamp("3Diter", 20,count,0);
1415 Float_t tempgridSizeR = gridSizeR * ione ;
1416 Float_t tempratioPhi = ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
1417 Float_t tempratioZ = ratioZ * ione * ione / ( jone * jone ) ;
1419 std::vector<float> coef1(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1420 std::vector<float> coef2(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1421 std::vector<float> coef3(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1422 std::vector<float> coef4(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
1424 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1425 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1426 coef1[i] = 1.0 + tempgridSizeR/(2*radius);
1427 coef2[i] = 1.0 - tempgridSizeR/(2*radius);
1428 coef3[i] = tempratioPhi/(radius*radius);
1429 coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
1432 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1433 TMatrixD &chargeDensity = *arrayofChargeDensities[m] ;
1434 TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
1435 for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
1436 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1437 for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
1438 if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
1439 else { // Add up all enclosed charge density contributions within 1/2 unit in all directions
1440 Float_t weight = 0.0 ;
1442 sumChargeDensity(i,j) = 0.0 ;
1443 for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
1444 for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
1445 if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
1448 sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
1449 sum += weight*radius ;
1452 sumChargeDensity(i,j) /= sum ;
1454 sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
1459 for ( Int_t k = 1 ; k <= iterations; k++ ) {
1461 // over-relaxation index, >= 1 but < 2
1462 Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
1463 Float_t overRelaxM1 = overRelax - 1.0 ;
1465 std::vector<float> overRelaxcoef4(rows) ; // Do this the standard C++ way to avoid gcc extensions
1466 std::vector<float> overRelaxcoef5(rows) ; // Do this the standard C++ way to avoid gcc extensions
1468 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1469 overRelaxcoef4[i] = overRelax * coef4[i] ;
1470 overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
1473 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1475 mplus = m + 1; signplus = 1 ;
1476 mminus = m - 1 ; signminus = 1 ;
1477 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1478 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1479 if ( mminus < 0 ) mminus = 1 ;
1481 else if (symmetry==-1) { // Anti-symmetry in phi
1482 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1483 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1485 else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1486 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1487 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1489 TMatrixD& arrayV = *arrayofArrayV[m] ;
1490 TMatrixD& arrayVP = *arrayofArrayV[mplus] ;
1491 TMatrixD& arrayVM = *arrayofArrayV[mminus] ;
1492 TMatrixD& sumChargeDensity = *arrayofSumChargeDensities[m] ;
1493 Double_t *arrayVfast = arrayV.GetMatrixArray();
1494 Double_t *arrayVPfast = arrayVP.GetMatrixArray();
1495 Double_t *arrayVMfast = arrayVM.GetMatrixArray();
1496 Double_t *sumChargeDensityFast=sumChargeDensity.GetMatrixArray();
1499 // slow implementation
1500 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1501 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1503 arrayV(i,j) = ( coef2[i] * arrayV(i-ione,j)
1504 + tempratioZ * ( arrayV(i,j-jone) + arrayV(i,j+jone) )
1505 - overRelaxcoef5[i] * arrayV(i,j)
1506 + coef1[i] * arrayV(i+ione,j)
1507 + coef3[i] * ( signplus*arrayVP(i,j) + signminus*arrayVM(i,j) )
1508 + sumChargeDensity(i,j)
1509 ) * overRelaxcoef4[i] ;
1510 // Note: over-relax the solution at each step. This speeds up the convergance.
1514 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1515 Double_t *arrayVfastI = &(arrayVfast[i*columns]);
1516 Double_t *arrayVPfastI = &(arrayVPfast[i*columns]);
1517 Double_t *arrayVMfastI = &(arrayVMfast[i*columns]);
1518 Double_t *sumChargeDensityFastI=&(sumChargeDensityFast[i*columns]);
1519 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1520 Double_t /*resSlow*/resFast;
1521 // resSlow = ( coef2[i] * arrayV(i-ione,j)
1522 // + tempratioZ * ( arrayV(i,j-jone) + arrayV(i,j+jone) )
1523 // - overRelaxcoef5[i] * arrayV(i,j)
1524 // + coef1[i] * arrayV(i+ione,j)
1525 // + coef3[i] * ( signplus*arrayVP(i,j) + signminus*arrayVM(i,j) )
1526 // + sumChargeDensity(i,j)
1527 // ) * overRelaxcoef4[i] ;
1528 resFast = ( coef2[i] * arrayVfastI[j-columns*ione]
1529 + tempratioZ * ( arrayVfastI[j-jone] + arrayVfastI[j+jone] )
1530 - overRelaxcoef5[i] * arrayVfastI[j]
1531 + coef1[i] * arrayVfastI[j+columns*ione]
1532 + coef3[i] * ( signplus* arrayVPfastI[j] + signminus*arrayVMfastI[j])
1533 + sumChargeDensityFastI[j]
1534 ) * overRelaxcoef4[i] ;
1535 // if (resSlow!=resFast){
1536 // printf("problem\t%d\t%d\t%f\t%f\t%f\n",i,j,resFast,resSlow,resFast-resSlow);
1538 arrayVfastI[j]=resFast;
1539 // Note: over-relax the solution at each step. This speeds up the convergance.
1544 if ( k == iterations ) { // After full solution is achieved, copy low resolution solution into higher res array
1545 for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
1546 for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
1549 arrayV(i+ione/2,j) = ( arrayV(i+ione,j) + arrayV(i,j) ) / 2 ;
1550 if ( i == ione ) arrayV(i-ione/2,j) = ( arrayV(0,j) + arrayV(ione,j) ) / 2 ;
1553 arrayV(i,j+jone/2) = ( arrayV(i,j+jone) + arrayV(i,j) ) / 2 ;
1554 if ( j == jone ) arrayV(i,j-jone/2) = ( arrayV(i,0) + arrayV(i,jone) ) / 2 ;
1556 if ( ione > 1 && jone > 1 ) {
1557 arrayV(i+ione/2,j+jone/2) = ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
1558 if ( i == ione ) arrayV(i-ione/2,j-jone/2) = ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
1559 if ( j == jone ) arrayV(i-ione/2,j-jone/2) = ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
1560 // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
1569 ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
1570 jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
1574 //Differentiate V(r) and solve for E(r) using special equations for the first and last row
1575 //Integrate E(r)/E(z) from point of origin to pad plane
1576 AliSysInfo::AddStamp("CalcField", 100,0,0);
1578 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1579 TMatrixD& arrayV = *arrayofArrayV[m] ;
1580 TMatrixD& eroverEz = *arrayofEroverEz[m] ;
1582 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1584 // Differentiate in R
1585 for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
1586 arrayE(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
1587 arrayE(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
1589 for ( Int_t i = 0 ; i < rows ; i++ ) {
1590 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1591 eroverEz(i,j) = 0.0 ;
1592 for ( Int_t k = j ; k < columns ; k++ ) {
1594 eroverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1595 if ( index != 4 ) index = 4; else index = 2 ;
1597 if ( index == 4 ) eroverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1598 if ( index == 2 ) eroverEz(i,j) +=
1599 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1600 if ( j == columns-2 ) eroverEz(i,j) =
1601 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1602 if ( j == columns-1 ) eroverEz(i,j) = 0.0 ;
1605 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1606 // eroverEz.Draw("surf") ; } // JT test
1608 AliSysInfo::AddStamp("IntegrateEr", 120,0,0);
1610 //Differentiate V(r) and solve for E(phi)
1611 //Integrate E(phi)/E(z) from point of origin to pad plane
1613 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1615 mplus = m + 1; signplus = 1 ;
1616 mminus = m - 1 ; signminus = 1 ;
1617 if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1618 if ( mplus > phislices-1 ) mplus = phislices - 2 ;
1619 if ( mminus < 0 ) mminus = 1 ;
1621 else if (symmetry==-1) { // Anti-symmetry in phi
1622 if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
1623 if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
1625 else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
1626 if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
1627 if ( mminus < 0 ) mminus = m - 1 + phislices ;
1629 TMatrixD &arrayVP = *arrayofArrayV[mplus] ;
1630 TMatrixD &arrayVM = *arrayofArrayV[mminus] ;
1631 TMatrixD &ePhioverEz = *arrayofEPhioverEz[m] ;
1632 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1633 // Differentiate in Phi
1634 for ( Int_t i = 0 ; i < rows ; i++ ) {
1635 Float_t radius = fgkIFCRadius + i*gridSizeR ;
1636 arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
1639 for ( Int_t i = 0 ; i < rows ; i++ ) {
1640 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1641 ePhioverEz(i,j) = 0.0 ;
1642 for ( Int_t k = j ; k < columns ; k++ ) {
1644 ePhioverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
1645 if ( index != 4 ) index = 4; else index = 2 ;
1647 if ( index == 4 ) ePhioverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
1648 if ( index == 2 ) ePhioverEz(i,j) +=
1649 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
1650 if ( j == columns-2 ) ePhioverEz(i,j) =
1651 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
1652 if ( j == columns-1 ) ePhioverEz(i,j) = 0.0 ;
1655 // if ( m == 5 ) { TCanvas* c2 = new TCanvas("arrayE","arrayE",50,50,840,600) ; c2 -> cd() ;
1656 // arrayE.Draw("surf") ; } // JT test
1658 AliSysInfo::AddStamp("IntegrateEphi", 130,0,0);
1661 // Differentiate V(r) and solve for E(z) using special equations for the first and last row
1662 // Integrate (E(z)-Ezstd) from point of origin to pad plane
1664 for ( Int_t m = 0 ; m < phislices ; m++ ) {
1665 TMatrixD& arrayV = *arrayofArrayV[m] ;
1666 TMatrixD& deltaEz = *arrayofDeltaEz[m] ;
1668 // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
1669 for ( Int_t i = 0 ; i < rows ; i++) {
1670 for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
1671 arrayE(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
1672 arrayE(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
1675 for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
1677 for ( Int_t i = 0 ; i < rows ; i++ ) {
1678 Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
1679 deltaEz(i,j) = 0.0 ;
1680 for ( Int_t k = j ; k < columns ; k++ ) {
1681 deltaEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k) ;
1682 if ( index != 4 ) index = 4; else index = 2 ;
1684 if ( index == 4 ) deltaEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1) ;
1685 if ( index == 2 ) deltaEz(i,j) +=
1686 (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
1687 if ( j == columns-2 ) deltaEz(i,j) =
1688 (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
1689 if ( j == columns-1 ) deltaEz(i,j) = 0.0 ;
1693 // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
1694 // eroverEz.Draw("surf") ; } // JT test
1696 // calculate z distortion from the integrated Delta Ez residuals
1697 // and include the aquivalence (Volt to cm) of the ROC shift !!
1699 for ( Int_t j = 0 ; j < columns ; j++ ) {
1700 for ( Int_t i = 0 ; i < rows ; i++ ) {
1702 // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
1703 deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
1705 // ROC Potential in cm aquivalent
1706 Double_t dzROCShift = arrayV(i, columns -1)/ezField;
1707 if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift; // add the ROC misaligment
1712 } // end loop over phi
1713 AliSysInfo::AddStamp("IntegrateEz", 140,0,0);
1716 for ( Int_t k = 0 ; k < phislices ; k++ )
1718 arrayofSumChargeDensities[k]->Delete() ;
1727 Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
1729 // Helperfunction: Check if integer is a power of 2
1732 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
1733 if ( j == 1 ) return(1) ; // True
1734 return(0) ; // False
1738 AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
1740 // Fit the track parameters - without and with distortion
1741 // 1. Space points in the TPC are simulated along the trajectory
1742 // 2. Space points distorted
1743 // 3. Fits the non distorted and distroted track to the reference plane at refX
1744 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
1746 // trackIn - input track parameters
1747 // refX - reference X to fit the track
1748 // dir - direction - out=1 or in=-1
1749 // pcstream - debug streamer to check the results
1751 // see AliExternalTrackParam.h documentation:
1752 // track1.fP[0] - local y (rphi)
1754 // track1.fP[2] - sinus of local inclination angle
1755 // track1.fP[3] - tangent of deep angle
1756 // track1.fP[4] - 1/pt
1758 AliTPCROC * roc = AliTPCROC::Instance();
1759 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
1760 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
1761 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1762 const Double_t kMaxSnp = 0.85;
1763 const Double_t kSigmaY=0.1;
1764 const Double_t kSigmaZ=0.1;
1765 const Double_t kMaxR=500;
1766 const Double_t kMaxZ=500;
1768 const Double_t kMaxZ0=220;
1769 const Double_t kZcut=3;
1770 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1774 AliExternalTrackParam track(trackIn); //
1776 AliTrackPointArray pointArray0(npoints0);
1777 AliTrackPointArray pointArray1(npoints0);
1779 if (!AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,5,kTRUE,kMaxSnp)) return 0;
1781 // simulate the track
1783 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
1784 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
1785 if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) return 0;
1787 xyz[0]+=gRandom->Gaus(0,0.000005);
1788 xyz[1]+=gRandom->Gaus(0,0.000005);
1789 xyz[2]+=gRandom->Gaus(0,0.000005);
1790 if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
1791 if (TMath::Abs(track.GetX())<kRTPC0) continue;
1792 if (TMath::Abs(track.GetX())>kRTPC1) continue;
1793 AliTrackPoint pIn0; // space point
1795 Int_t sector= (xyz[2]>0)? 0:18;
1796 pointArray0.GetPoint(pIn0,npoints);
1797 pointArray1.GetPoint(pIn1,npoints);
1798 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
1799 Float_t distPoint[3]={static_cast<Float_t>(xyz[0]),static_cast<Float_t>(xyz[1]),static_cast<Float_t>(xyz[2])};
1800 DistortPoint(distPoint, sector);
1801 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
1802 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
1804 track.Rotate(alpha);
1805 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
1806 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
1807 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
1808 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
1809 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
1810 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
1811 pointArray0.AddPoint(npoints, &pIn0);
1812 pointArray1.AddPoint(npoints, &pIn1);
1814 if (npoints>=npoints0) break;
1816 if (npoints<npoints0/4.) return 0;
1820 AliExternalTrackParam *track0=0;
1821 AliExternalTrackParam *track1=0;
1822 AliTrackPoint point1,point2,point3;
1823 if (dir==1) { //make seed inner
1824 pointArray0.GetPoint(point1,1);
1825 pointArray0.GetPoint(point2,11);
1826 pointArray0.GetPoint(point3,21);
1828 if (dir==-1){ //make seed outer
1829 pointArray0.GetPoint(point1,npoints-21);
1830 pointArray0.GetPoint(point2,npoints-11);
1831 pointArray0.GetPoint(point3,npoints-1);
1833 if ((TMath::Abs(point1.GetX()-point3.GetX())+TMath::Abs(point1.GetY()-point3.GetY()))<10){
1834 printf("fit points not properly initialized\n");
1837 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
1838 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
1839 track0->ResetCovariance(10);
1840 track1->ResetCovariance(10);
1841 if (TMath::Abs(AliTrackerBase::GetBz())<0.01){
1842 ((Double_t*)track0->GetParameter())[4]= trackIn.GetParameter()[4];
1843 ((Double_t*)track1->GetParameter())[4]= trackIn.GetParameter()[4];
1845 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
1846 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
1850 pointArray0.GetPoint(pIn0,ipoint);
1851 pointArray1.GetPoint(pIn1,ipoint);
1852 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
1853 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
1854 if (TMath::Abs(prot0.GetX())<kRTPC0) continue;
1855 if (TMath::Abs(prot0.GetX())>kRTPC1) continue;
1857 if (!AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
1858 if (!AliTrackerBase::PropagateTrackTo(track1,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
1859 if (TMath::Abs(track0->GetZ())>kMaxZ) break;
1860 if (TMath::Abs(track0->GetX())>kMaxR) break;
1861 if (TMath::Abs(track1->GetZ())>kMaxZ) break;
1862 if (TMath::Abs(track1->GetX())>kMaxR) break;
1863 if (dir>0 && track1->GetX()>refX) continue;
1864 if (dir<0 && track1->GetX()<refX) continue;
1865 if (TMath::Abs(track1->GetZ())<kZcut)continue;
1866 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
1868 Double_t pointPos[2]={0,0};
1869 Double_t pointCov[3]={0,0,0};
1870 pointPos[0]=prot0.GetY();//local y
1871 pointPos[1]=prot0.GetZ();//local z
1872 pointCov[0]=prot0.GetCov()[3];//simay^2
1873 pointCov[1]=prot0.GetCov()[4];//sigmayz
1874 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
1875 if (!track0->Update(pointPos,pointCov)) break;
1877 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
1878 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
1879 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
1881 pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
1882 pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
1883 pointCov[0]=prot1.GetCov()[3];//simay^2
1884 pointCov[1]=prot1.GetCov()[4];//sigmayz
1885 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
1886 if (!track1->Update(pointPos,pointCov)) break;
1890 if (npoints2<npoints/4.) return 0;
1891 AliTrackerBase::PropagateTrackTo(track0,refX,kMass,5.,kTRUE,kMaxSnp);
1892 AliTrackerBase::PropagateTrackTo(track0,refX,kMass,1.,kTRUE,kMaxSnp);
1893 track1->Rotate(track0->GetAlpha());
1894 AliTrackerBase::PropagateTrackTo(track1,track0->GetX(),kMass,5.,kFALSE,kMaxSnp);
1896 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
1897 "point0.="<<&pointArray0<< // points
1898 "point1.="<<&pointArray1<< // distorted points
1899 "trackIn.="<<&track<< // original track
1900 "track0.="<<track0<< // fitted track
1901 "track1.="<<track1<< // fitted distorted track
1903 new(&trackIn) AliExternalTrackParam(*track0);
1912 TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
1914 // create the distortion tree on a mesh with granularity given by step
1915 // return the tree with distortions at given position
1916 // Map is created on the mesh with given step size
1918 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
1920 for (Double_t x= -250; x<250; x+=step){
1921 for (Double_t y= -250; y<250; y+=step){
1922 Double_t r = TMath::Sqrt(x*x+y*y);
1924 if (r>250) continue;
1925 for (Double_t z= -250; z<250; z+=step){
1926 Int_t roc=(z>0)?0:18;
1930 Double_t phi = TMath::ATan2(y,x);
1931 DistortPoint(xyz,roc);
1932 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1933 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
1934 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
1935 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
1936 Double_t dx = xyz[0]-x;
1937 Double_t dy = xyz[1]-y;
1938 Double_t dz = xyz[2]-z;
1940 Double_t drphi=(phi1-phi)*r;
1941 (*pcstream)<<"distortion"<<
1942 "x="<<x<< // original position
1947 "x1="<<xyz[0]<< // distorted position
1953 "dx="<<dx<< // delta position
1963 TFile f(Form("correction%s.root",GetName()));
1964 TTree * tree = (TTree*)f.Get("distortion");
1965 TTree * tree2= tree->CopyTree("1");
1966 tree2->SetName(Form("dist%s",GetName()));
1967 tree2->SetDirectory(0);
1975 void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1978 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1979 // calculates partial distortions
1980 // Partial distortion is stored in the resulting tree
1981 // Output is storred in the file distortion_<dettype>_<partype>.root
1982 // Partial distortion is stored with the name given by correction name
1985 // Parameters of function:
1986 // input - input tree
1987 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex , 3 - TPC-TOF, 4 - TPCTPC track crossing
1988 // ppype - parameter type
1989 // corrArray - array with partial corrections
1990 // step - skipe entries - if 1 all entries processed - it is slow
1991 // debug 0 if debug on also space points dumped - it is slow
1993 const Double_t kMaxSnp = 0.85;
1994 const Double_t kcutSnp=0.25;
1995 const Double_t kcutTheta=1.;
1996 const Double_t kRadiusTPC=85;
1997 // AliTPCROC *tpcRoc =AliTPCROC::Instance();
1999 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2000 // const Double_t kB2C=-0.299792458e-3;
2001 const Int_t kMinEntries=20;
2002 Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
2005 tinput->SetBranchAddress("run",&run);
2006 tinput->SetBranchAddress("theta",&theta);
2007 tinput->SetBranchAddress("phi", &phi);
2008 tinput->SetBranchAddress("snp",&snp);
2009 tinput->SetBranchAddress("mean",&mean);
2010 tinput->SetBranchAddress("rms",&rms);
2011 tinput->SetBranchAddress("entries",&entries);
2012 tinput->SetBranchAddress("sector",§or);
2013 tinput->SetBranchAddress("dsec",&dsec);
2014 tinput->SetBranchAddress("refX",&refX);
2015 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
2017 Int_t nentries=tinput->GetEntries();
2018 Int_t ncorr=corrArray->GetEntries();
2019 Double_t corrections[100]={0}; //
2021 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
2023 if (dtype==5 || dtype==6) dtype=4;
2024 if (dtype==0) { dir=-1;}
2025 if (dtype==1) { dir=1;}
2026 if (dtype==2) { dir=-1;}
2027 if (dtype==3) { dir=1;}
2028 if (dtype==4) { dir=-1;}
2030 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
2031 tinput->GetEntry(ientry);
2032 if (TMath::Abs(snp)>kMaxSnp) continue;
2035 if (dtype==2) tPar[1]=theta*kRadiusTPC;
2038 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
2040 // tracks crossing CE
2041 tPar[1]=0; // track at the CE
2042 //if (TMath::Abs(theta) <0.05) continue; // deep cross
2045 if (TMath::Abs(snp) >kcutSnp) continue;
2046 if (TMath::Abs(theta) >kcutTheta) continue;
2047 printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
2048 Double_t bz=AliTrackerBase::GetBz();
2049 if (dtype !=4) { //exclude TPC - for TPC mainly non primary tracks
2050 if (dtype!=2 && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
2052 if (dtype==2 && TMath::Abs(bz)>0.1 ) {
2053 tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);//
2054 // snp at the TPC inner radius in case the vertex match used
2058 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
2059 AliExternalTrackParam track(refX,phi,tPar,cov);
2063 Double_t pt=1./tPar[4];
2064 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
2065 //if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used
2066 Double_t refXD=refX;
2067 (*pcstream)<<"fit"<<
2068 "run="<<run<< // run number
2069 "bz="<<bz<< // magnetic filed used
2070 "dtype="<<dtype<< // detector match type
2071 "ptype="<<ptype<< // parameter type
2072 "theta="<<theta<< // theta
2073 "phi="<<phi<< // phi
2074 "snp="<<snp<< // snp
2075 "mean="<<mean<< // mean dist value
2076 "rms="<<rms<< // rms
2079 "refX="<<refXD<< // referece X as double
2080 "gx="<<xyz[0]<< // global position at reference
2081 "gy="<<xyz[1]<< // global position at reference
2082 "gz="<<xyz[2]<< // global position at reference
2083 "dRrec="<<dRrec<< // delta Radius in reconstruction
2085 "id="<<id<< // track id
2086 "entries="<<entries;// number of entries in bin
2089 if (entries<kMinEntries) isOK=kFALSE;
2091 if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
2092 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
2093 corrections[icorr]=0;
2094 if (entries>kMinEntries){
2095 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
2096 AliExternalTrackParam *trackOut = 0;
2097 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
2098 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
2099 if (dtype==0) {dir= -1;}
2100 if (dtype==1) {dir= 1;}
2101 if (dtype==2) {dir= -1;}
2102 if (dtype==3) {dir= 1;}
2105 if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
2106 if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE;
2107 if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
2108 // trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
2110 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
2113 corrections[icorr]=0;
2116 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out
2118 (*pcstream)<<"fit"<<
2119 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
2122 if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
2124 // special case of the TPC tracks crossing the CE
2126 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
2127 corrections[icorr]=0;
2128 if (entries>kMinEntries){
2129 AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
2130 AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet
2131 AliExternalTrackParam *trackOut0 = 0;
2132 AliExternalTrackParam *trackOut1 = 0;
2134 if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
2135 if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
2136 if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
2137 if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
2139 if (trackOut0 && trackOut1){
2140 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
2141 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2142 if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2143 if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
2145 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
2146 if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2147 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2148 if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE;
2149 if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
2151 corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
2152 corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
2154 if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
2155 (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
2156 (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
2157 (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
2158 (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
2159 (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
2166 corrections[icorr]=0;
2170 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out no in lookup
2172 (*pcstream)<<"fit"<<
2173 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
2176 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
2185 void AliTPCCorrection::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
2188 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
2189 // calculates partial distortions
2190 // Partial distortion is stored in the resulting tree
2191 // Output is storred in the file distortion_<dettype>_<partype>.root
2192 // Partial distortion is stored with the name given by correction name
2195 // Parameters of function:
2196 // input - input tree
2197 // dtype - distortion type 10 - IROC-OROC
2198 // ppype - parameter type
2199 // corrArray - array with partial corrections
2200 // step - skipe entries - if 1 all entries processed - it is slow
2201 // debug 0 if debug on also space points dumped - it is slow
2203 const Double_t kMaxSnp = 0.8;
2204 const Int_t kMinEntries=200;
2205 // AliTPCROC *tpcRoc =AliTPCROC::Instance();
2207 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2208 // const Double_t kB2C=-0.299792458e-3;
2209 Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
2214 tinput->SetBranchAddress("run",&run);
2215 tinput->SetBranchAddress("theta",&theta);
2216 tinput->SetBranchAddress("phi", &phi);
2217 tinput->SetBranchAddress("snp",&snp);
2218 tinput->SetBranchAddress("mean",&mean);
2219 tinput->SetBranchAddress("rms",&rms);
2220 tinput->SetBranchAddress("entries",&entries);
2221 tinput->SetBranchAddress("sector",§or);
2222 tinput->SetBranchAddress("dsec",&dsec);
2223 tinput->SetBranchAddress("refX",&refXD);
2224 tinput->SetBranchAddress("z",&globalZ);
2225 tinput->SetBranchAddress("isec0",&isec0);
2226 tinput->SetBranchAddress("isec1",&isec1);
2227 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
2229 Int_t nentries=tinput->GetEntries();
2230 Int_t ncorr=corrArray->GetEntries();
2231 Double_t corrections[100]={0}; //
2233 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
2236 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
2237 tinput->GetEntry(ientry);
2240 if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1; // IROC-OROC - opposite side
2241 if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2; // IROC-OROC - same side
2242 if (dtype==10 && id==-1) continue;
2249 tPar[4]=(gRandom->Rndm()-0.1)*0.2; //
2250 Double_t pt=1./tPar[4];
2252 printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
2253 Double_t bz=AliTrackerBase::GetBz();
2254 AliExternalTrackParam track(refX,phi,tPar,cov);
2255 Double_t xyz[3],xyzIn[3],xyzOut[3];
2257 track.GetXYZAt(85,bz,xyzIn);
2258 track.GetXYZAt(245,bz,xyzOut);
2259 Double_t phiIn = TMath::ATan2(xyzIn[1],xyzIn[0]);
2260 Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
2261 Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
2262 Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
2263 Int_t sectorIn = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
2264 Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
2267 if (sectorIn!=sectorOut) isOK=kFALSE; // requironment - cluster in the same sector
2268 if (sectorIn!=sectorRef) isOK=kFALSE; // requironment - cluster in the same sector
2269 if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE; // requironment - minimal amount of tracks in bin
2271 if (TMath::Abs(theta)>1) isOK=kFALSE;
2273 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
2275 (*pcstream)<<"fit"<<
2277 "bz="<<bz<< // magnetic filed used
2278 "dtype="<<dtype<< // detector match type
2279 "ptype="<<ptype<< // parameter type
2280 "theta="<<theta<< // theta
2281 "phi="<<phi<< // phi
2282 "snp="<<snp<< // snp
2283 "mean="<<mean<< // mean dist value
2284 "rms="<<rms<< // rms
2287 "refX="<<refXD<< // referece X
2288 "gx="<<xyz[0]<< // global position at reference
2289 "gy="<<xyz[1]<< // global position at reference
2290 "gz="<<xyz[2]<< // global position at reference
2291 "dRrec="<<dRrec<< // delta Radius in reconstruction
2293 "id="<<id<< // track id
2294 "entries="<<entries;// number of entries in bin
2296 AliExternalTrackParam *trackOut0 = 0;
2297 AliExternalTrackParam *trackOut1 = 0;
2298 AliExternalTrackParam *ptrackIn0 = 0;
2299 AliExternalTrackParam *ptrackIn1 = 0;
2301 for (Int_t icorr=0; icorr<ncorr; icorr++) {
2303 // special case of the TPC tracks crossing the CE
2305 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
2306 corrections[icorr]=0;
2307 if (entries>kMinEntries &&isOK){
2308 AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
2309 AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
2310 ptrackIn1=&trackIn0;
2311 ptrackIn0=&trackIn1;
2313 if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
2314 if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
2315 if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
2316 if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
2318 if (trackOut0 && trackOut1){
2320 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp)) isOK=kFALSE;
2321 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2322 // rotate all tracks to the same frame
2323 if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2324 if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2325 if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
2327 if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2328 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2329 if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
2331 corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
2332 corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
2333 (*pcstream)<<"fitDebug"<< // just to debug the correction
2335 "pIn0.="<<ptrackIn0<<
2336 "pIn1.="<<ptrackIn1<<
2337 "pOut0.="<<trackOut0<<
2338 "pOut1.="<<trackOut1<<
2344 corrections[icorr]=0;
2348 (*pcstream)<<"fit"<<
2349 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
2352 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
2359 void AliTPCCorrection::MakeLaserDistortionTreeOld(TTree* tree, TObjArray *corrArray, Int_t itype){
2361 // Make a laser fit tree for global minimization
2363 const Double_t cutErrY=0.1;
2364 const Double_t cutErrZ=0.1;
2365 const Double_t kEpsilon=0.00000001;
2366 const Double_t kMaxDist=1.; // max distance - space correction
2367 const Double_t kMaxRMS=0.05; // max distance -between point and local mean
2372 AliTPCLaserTrack *ltr=0;
2373 AliTPCLaserTrack::LoadTracks();
2374 tree->SetBranchAddress("dY.",&vecdY);
2375 tree->SetBranchAddress("dZ.",&vecdZ);
2376 tree->SetBranchAddress("eY.",&veceY);
2377 tree->SetBranchAddress("eZ.",&veceZ);
2378 tree->SetBranchAddress("LTr.",<r);
2379 Int_t entries= tree->GetEntries();
2380 TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
2381 Double_t bz=AliTrackerBase::GetBz();
2384 for (Int_t ientry=0; ientry<entries; ientry++){
2385 tree->GetEntry(ientry);
2386 if (!ltr->GetVecGX()){
2387 ltr->UpdatePoints();
2389 TVectorD * delta= (itype==0)? vecdY:vecdZ;
2390 TVectorD * err= (itype==0)? veceY:veceZ;
2391 TLinearFitter fitter(2,"pol1");
2392 for (Int_t iter=0; iter<2; iter++){
2393 Double_t kfit0=0, kfit1=0;
2394 Int_t npoints=fitter.GetNpoints();
2397 kfit0=fitter.GetParameter(0);
2398 kfit1=fitter.GetParameter(1);
2400 for (Int_t irow=0; irow<159; irow++){
2403 Int_t nentries = 1000;
2404 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
2405 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
2408 Int_t first3=TMath::Max(irow-3,0);
2409 Int_t last3 =TMath::Min(irow+3,159);
2411 if ((*ltr->GetVecSec())[irow]>=0 && err) {
2412 for (Int_t jrow=first3; jrow<=last3; jrow++){
2413 if ((*ltr->GetVecSec())[irow]!= (*ltr->GetVecSec())[jrow]) continue;
2414 if ((*err)[jrow]<kEpsilon) continue;
2415 array[counter]=(*delta)[jrow];
2422 rms3 = TMath::RMS(counter,array);
2423 mean3 = TMath::Mean(counter,array);
2427 Double_t phi =(*ltr->GetVecPhi())[irow];
2428 Double_t theta =ltr->GetTgl();
2429 Double_t mean=delta->GetMatrixArray()[irow];
2430 Double_t gx=0,gy=0,gz=0;
2431 Double_t snp = (*ltr->GetVecP2())[irow];
2433 // Double_t rms = err->GetMatrixArray()[irow];
2435 gx = (*ltr->GetVecGX())[irow];
2436 gy = (*ltr->GetVecGY())[irow];
2437 gz = (*ltr->GetVecGZ())[irow];
2439 // get delta R used in reconstruction
2440 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
2441 AliTPCCorrection * correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
2442 // const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
2443 //Double_t xyz0[3]={gx,gy,gz};
2444 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
2445 Double_t fphi = TMath::ATan2(gy,gx);
2446 Double_t fsector = 9.*fphi/TMath::Pi();
2447 if (fsector<0) fsector+=18;
2448 Double_t dsec = fsector-Int_t(fsector)-0.5;
2450 Int_t id= ltr->GetId();
2454 Float_t xyz1[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
2455 Int_t sector=(gz>0)?0:18;
2456 correction->CorrectPoint(xyz1, sector);
2457 refX=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
2460 if (TMath::Abs(rms3)>kMaxRMS) isOK=kFALSE;
2461 if (TMath::Abs(mean-mean3)>kMaxRMS) isOK=kFALSE;
2462 if (counter<4) isOK=kFALSE;
2463 if (npoints<90) isOK=kFALSE;
2465 fitter.AddPoint(&refX,mean);
2467 Double_t deltaF=kfit0+kfit1*refX;
2469 (*pcstream)<<"fitFull"<< // dumpe also intermediate results
2470 "bz="<<bz<< // magnetic filed used
2471 "dtype="<<dtype<< // detector match type
2472 "ptype="<<itype<< // parameter type
2473 "theta="<<theta<< // theta
2474 "phi="<<phi<< // phi
2475 "snp="<<snp<< // snp
2476 "mean="<<mean3<< // mean dist value
2477 "rms="<<rms3<< // rms
2479 "npoints="<<npoints<< //number of points
2480 "mean3="<<mean3<< // mean dist value
2481 "rms3="<<rms3<< // rms
2482 "counter="<<counter<<
2483 "sector="<<fsector<<
2486 "refX="<<refX<< // reference radius
2487 "gx="<<gx<< // global position
2488 "gy="<<gy<< // global position
2489 "gz="<<gz<< // global position
2490 "dRrec="<<dRrec<< // delta Radius in reconstruction
2491 "id="<<id<< //bundle
2492 "entries="<<nentries<<// number of entries in bin
2495 if (iter==1) (*pcstream)<<"fit"<< // dump valus for fit
2496 "bz="<<bz<< // magnetic filed used
2497 "dtype="<<dtype<< // detector match type
2498 "ptype="<<itype<< // parameter type
2499 "theta="<<theta<< // theta
2500 "phi="<<phi<< // phi
2501 "snp="<<snp<< // snp
2502 "mean="<<mean3<< // mean dist value
2503 "rms="<<rms3<< // rms
2504 "sector="<<fsector<<
2507 "refX="<<refX<< // reference radius
2508 "gx="<<gx<< // global position
2509 "gy="<<gy<< // global position
2510 "gz="<<gz<< // global position
2511 "dRrec="<<dRrec<< // delta Radius in reconstruction
2513 "id="<<id<< //bundle
2514 "entries="<<nentries;// number of entries in bin
2517 Double_t ky = TMath::Tan(TMath::ASin(snp));
2518 Int_t ncorr = corrArray->GetEntries();
2519 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
2520 Double_t phi0 = TMath::ATan2(gy,gx);
2521 Double_t distortions[1000]={0};
2522 Double_t distortionsR[1000]={0};
2524 for (Int_t icorr=0; icorr<ncorr; icorr++) {
2525 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
2526 Float_t distPoint[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
2527 Int_t sector= (gz>0)? 0:18;
2529 corr->DistortPoint(distPoint, sector);
2531 // Double_t value=distPoint[2]-gz;
2532 if (itype==0 && r0>1){
2533 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
2534 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
2535 Double_t drphi= r0*(phi1-phi0);
2536 Double_t dr = r1-r0;
2537 distortions[icorr] = drphi-ky*dr;
2538 distortionsR[icorr] = dr;
2540 if (TMath::Abs(distortions[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE; }
2541 if (TMath::Abs(distortionsR[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE;}
2542 (*pcstream)<<"fit"<<
2543 Form("%s=",corr->GetName())<<distortions[icorr]; // dump correction value
2545 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
2555 void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type, Int_t integ){
2557 // make a distortion map out ou fthe residual histogram
2558 // Results are written to the debug streamer - pcstream
2560 // his0 - input (4D) residual histogram
2561 // pcstream - file to write the tree
2563 // refX - track matching reference X
2564 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
2566 // OBJ: TAxis #Delta #Delta
2567 // OBJ: TAxis tanTheta tan(#Theta)
2568 // OBJ: TAxis phi #phi
2569 // OBJ: TAxis snp snp
2571 // marian.ivanov@cern.ch
2572 const Int_t kMinEntries=10;
2573 Double_t bz=AliTrackerBase::GetBz();
2574 Int_t idim[4]={0,1,2,3};
2578 Int_t nbins3=his0->GetAxis(3)->GetNbins();
2579 Int_t first3=his0->GetAxis(3)->GetFirst();
2580 Int_t last3 =his0->GetAxis(3)->GetLast();
2582 for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - local angle
2583 his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3));
2584 Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3);
2585 THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3
2587 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
2588 Int_t first2 = his3->GetAxis(2)->GetFirst();
2589 Int_t last2 = his3->GetAxis(2)->GetLast();
2591 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi
2592 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2));
2593 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
2594 THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
2595 Int_t nbins1 = his2->GetAxis(1)->GetNbins();
2596 Int_t first1 = his2->GetAxis(1)->GetFirst();
2597 Int_t last1 = his2->GetAxis(1)->GetLast();
2598 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
2600 Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1);
2601 his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
2602 if (TMath::Abs(x1)<0.1){
2603 if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1));
2604 if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1));
2606 if (TMath::Abs(x1)<0.06){
2607 his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
2609 TH1 * hisDelta = his2->Projection(0);
2611 Double_t entries = hisDelta->GetEntries();
2612 Double_t mean=0, rms=0;
2613 if (entries>kMinEntries){
2614 mean = hisDelta->GetMean();
2615 rms = hisDelta->GetRMS();
2617 Double_t sector = 9.*x2/TMath::Pi();
2618 if (sector<0) sector+=18;
2619 Double_t dsec = sector-Int_t(sector)-0.5;
2621 (*pcstream)<<hname<<
2626 "z="<<z<< // dummy z
2628 "entries="<<entries<<
2631 "refX="<<refX<< // track matching refernce plane
2637 //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
2648 void AliTPCCorrection::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){
2650 // make a distortion map out ou fthe residual histogram
2651 // Results are written to the debug streamer - pcstream
2653 // his0 - input (4D) residual histogram
2654 // pcstream - file to write the tree
2656 // refX - track matching reference X
2657 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
2658 // marian.ivanov@cern.ch
2661 // Collection name='TObjArray', class='TObjArray', size=16
2662 // 0. OBJ: TAxis #Delta #Delta
2663 // 1. OBJ: TAxis N_{cl} N_{cl}
2664 // 2. OBJ: TAxis dca_{r} (cm) dca_{r} (cm)
2665 // 3. OBJ: TAxis z (cm) z (cm)
2666 // 4. OBJ: TAxis sin(#phi) sin(#phi)
2667 // 5. OBJ: TAxis tan(#theta) tan(#theta)
2668 // 6. OBJ: TAxis 1/pt (1/GeV) 1/pt (1/GeV)
2669 // 7. OBJ: TAxis pt (GeV) pt (GeV)
2670 // 8. OBJ: TAxis alpha alpha
2671 const Int_t kMinEntries=10;
2673 // 1. make default selections
2676 Int_t idim0[4]={0 , 5, 8, 3}; // delta, theta, alpha, z
2677 hisInput->GetAxis(1)->SetRangeUser(110,190); //long tracks
2678 hisInput->GetAxis(2)->SetRangeUser(-10,35); //tracks close to beam pipe
2679 hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance
2680 hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks"
2681 hisDelta= hisInput->Projection(0);
2682 hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS());
2684 THnSparse *his0= hisInput->Projection(4,idim0);
2686 // 2. Get mean in diferent bins
2688 Int_t nbins1=his0->GetAxis(1)->GetNbins();
2689 Int_t first1=his0->GetAxis(1)->GetFirst();
2690 Int_t last1 =his0->GetAxis(1)->GetLast();
2692 Double_t bz=AliTrackerBase::GetBz();
2693 Int_t idim[4]={0,1, 2, 3}; // delta, theta,alpha,z
2695 for (Int_t ibin1=first1; ibin1<=last1; ibin1++){ //axis 1 - theta
2697 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
2698 his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
2700 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
2701 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
2702 Int_t first3 = his1->GetAxis(3)->GetFirst();
2703 Int_t last3 = his1->GetAxis(3)->GetLast();
2705 for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - z at "vertex"
2706 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
2707 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
2709 his1->GetAxis(3)->SetRangeUser(-1,1);
2712 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
2713 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
2714 Int_t first2 = his3->GetAxis(2)->GetFirst();
2715 Int_t last2 = his3->GetAxis(2)->GetLast();
2717 for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){
2718 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
2719 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
2720 hisDelta = his3->Projection(0);
2722 Double_t entries = hisDelta->GetEntries();
2723 Double_t mean=0, rms=0;
2724 if (entries>kMinEntries){
2725 mean = hisDelta->GetMean();
2726 rms = hisDelta->GetRMS();
2728 Double_t sector = 9.*x2/TMath::Pi();
2729 if (sector<0) sector+=18;
2730 Double_t dsec = sector-Int_t(sector)-0.5;
2731 Double_t snp=0; // dummy snp - equal 0
2732 (*pcstream)<<hname<<
2734 "bz="<<bz<< // magnetic field
2735 "theta="<<x1<< // theta
2736 "phi="<<x2<< // phi (alpha)
2737 "z="<<x3<< // z at "vertex"
2738 "snp="<<snp<< // dummy snp
2739 "entries="<<entries<< // entries in bin
2740 "mean="<<mean<< // mean
2742 "refX="<<refX<< // track matching refernce plane
2743 "type="<<type<< // parameter type
2744 "sector="<<sector<< // sector
2745 "dsec="<<dsec<< // dummy delta sector
2748 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
2759 void AliTPCCorrection::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type){
2761 // make a distortion map out of the residual histogram
2762 // Results are written to the debug streamer - pcstream
2764 // his0 - input (4D) residual histogram
2765 // pcstream - file to write the tree
2767 // type - 0- y 1-z,2 -snp, 3-theta
2768 // marian.ivanov@cern.ch
2770 //Collection name='TObjArray', class='TObjArray', size=16
2771 //0 OBJ: TAxis delta delta
2772 //1 OBJ: TAxis phi phi
2773 //2 OBJ: TAxis localX localX
2774 //3 OBJ: TAxis kY kY
2775 //4 OBJ: TAxis kZ kZ
2776 //5 OBJ: TAxis is1 is1
2777 //6 OBJ: TAxis is0 is0
2779 //8. OBJ: TAxis IsPrimary IsPrimary
2781 const Int_t kMinEntries=10;
2782 THnSparse * hisSector0=0;
2783 TH1 * htemp=0; // histogram to calculate mean value of parameter
2784 Double_t bz=AliTrackerBase::GetBz();
2787 // Loop over pair of sector:
2798 for (Int_t isec0=0; isec0<72; isec0++){
2799 Int_t index0[9]={0, 4, 3, 7, 1, 2, 5, 6,8}; //regroup indeces
2801 //hisInput->GetAxis(8)->SetRangeUser(-0.1,0.4); // select secondaries only ? - get out later ?
2802 hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1);
2803 hisSector0=hisInput->Projection(7,index0);
2806 for (Int_t isec1=isec0+1; isec1<72; isec1++){
2807 //if (isec1!=isec0+36) continue;
2808 if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue;
2809 printf("Sectors %d\t%d\n",isec1,isec0);
2810 hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1);
2811 TH1 * hisX=hisSector0->Projection(5);
2812 Double_t refX= hisX->GetMean();
2814 TH1 *hisDelta=hisSector0->Projection(0);
2815 Double_t dmean = hisDelta->GetMean();
2816 Double_t drms = hisDelta->GetRMS();
2817 hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms);
2820 // 1. make default selections
2822 Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi }
2823 THnSparse *hisSector1= hisSector0->Projection(5,idim0);
2825 // 2. Get mean in diferent bins
2827 Int_t idim[5]={0, 1, 2, 3, 4}; // {delta, theta-1,snp-2 ,z-3, phi-4}
2829 // Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins();
2830 Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst();
2831 Int_t lastPhi =hisSector1->GetAxis(4)->GetLast();
2833 for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=1){ //axis 4 - phi
2837 Double_t xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi);
2838 Double_t psec = (9*xPhi/TMath::Pi());
2839 if (psec<0) psec+=18;
2840 Bool_t isOK0=kFALSE;
2841 Bool_t isOK1=kFALSE;
2842 if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.) isOK0=kTRUE;
2843 if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.) isOK1=kTRUE;
2844 if (!isOK0) continue;
2845 if (!isOK1) continue;
2847 hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi));
2848 if (isec1!=isec0+36) {
2849 hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi));
2852 htemp = hisSector1->Projection(4);
2853 xPhi=htemp->GetMean();
2855 THnSparse * hisPhi = hisSector1->Projection(4,idim);
2856 //Int_t nbinsZ = hisPhi->GetAxis(3)->GetNbins();
2857 Int_t firstZ = hisPhi->GetAxis(3)->GetFirst();
2858 Int_t lastZ = hisPhi->GetAxis(3)->GetLast();
2860 for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=1){ // axis 3 - z
2864 hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ));
2865 if (isec1!=isec0+36) {
2866 hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ));
2868 htemp = hisPhi->Projection(3);
2869 Double_t xZ= htemp->GetMean();
2871 THnSparse * hisZ= hisPhi->Projection(3,idim);
2872 //projected histogram according selection 3 -z
2875 //Int_t nbinsSnp = hisZ->GetAxis(2)->GetNbins();
2876 Int_t firstSnp = hisZ->GetAxis(2)->GetFirst();
2877 Int_t lastSnp = hisZ->GetAxis(2)->GetLast();
2878 for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){ // axis 2 - snp
2882 hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp));
2883 if (isec1!=isec0+36) {
2884 hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp));
2886 htemp = hisZ->Projection(2);
2887 Double_t xSnp= htemp->GetMean();
2889 THnSparse * hisSnp= hisZ->Projection(2,idim);
2890 //projected histogram according selection 2 - snp
2892 //Int_t nbinsTheta = hisSnp->GetAxis(1)->GetNbins();
2893 Int_t firstTheta = hisSnp->GetAxis(1)->GetFirst();
2894 Int_t lastTheta = hisSnp->GetAxis(1)->GetLast();
2896 for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){ // axis1 theta
2899 hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta));
2900 if (isec1!=isec0+36) {
2901 hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta));
2903 htemp = hisSnp->Projection(1);
2904 Double_t xTheta=htemp->GetMean();
2906 hisDelta = hisSnp->Projection(0);
2908 Double_t entries = hisDelta->GetEntries();
2909 Double_t mean=0, rms=0;
2910 if (entries>kMinEntries){
2911 mean = hisDelta->GetMean();
2912 rms = hisDelta->GetRMS();
2914 Double_t sector = 9.*xPhi/TMath::Pi();
2915 if (sector<0) sector+=18;
2916 Double_t dsec = sector-Int_t(sector)-0.5;
2917 Int_t dtype=1; // TPC alignment type
2918 (*pcstream)<<hname<<
2920 "bz="<<bz<< // magnetic field
2921 "ptype="<<type<< // parameter type
2922 "dtype="<<dtype<< // parameter type
2923 "isec0="<<isec0<< // sector 0
2924 "isec1="<<isec1<< // sector 1
2925 "sector="<<sector<< // sector as float
2926 "dsec="<<dsec<< // delta sector
2928 "theta="<<xTheta<< // theta
2929 "phi="<<xPhi<< // phi (alpha)
2931 "snp="<<xSnp<< // snp
2933 "entries="<<entries<< // entries in bin
2934 "mean="<<mean<< // mean
2935 "rms="<<rms<< // rms
2936 "refX="<<refX<< // track matching reference plane
2939 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);
2960 void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
2962 // Store object in the OCDB
2963 // By default the object is stored in the current directory
2964 // default comment consit of user name and the date
2966 TString ocdbStorage="";
2967 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
2968 AliCDBMetaData *metaData= new AliCDBMetaData();
2969 metaData->SetObjectClassName("AliTPCCorrection");
2970 metaData->SetResponsible("Marian Ivanov");
2971 metaData->SetBeamPeriod(1);
2972 metaData->SetAliRootVersion("05-25-01"); //root version
2973 TString userName=gSystem->GetFromPipe("echo $USER");
2974 TString date=gSystem->GetFromPipe("date");
2976 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
2977 if (comment) metaData->SetComment(comment);
2979 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
2980 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
2981 gStorage->Put(this, (*id1), metaData);
2985 void AliTPCCorrection::FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts){
2987 // Fast method to simulate the influence of the given distortion on the vertex reconstruction
2990 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
2991 if (!magF) AliError("Magneticd field - not initialized");
2992 Double_t bz = magF->SolenoidField(); //field in kGauss
2993 printf("bz: %f\n",bz);
2994 AliVertexerTracks *vertexer = new AliVertexerTracks(bz); // bz in kGauss
2996 TObjArray aTrk; // Original Track array of Aside
2997 TObjArray daTrk; // Distorted Track array of A side
2998 UShort_t *aId = new UShort_t[nTracks]; // A side Track ID
3001 UShort_t *cId = new UShort_t [nTracks];
3003 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
3004 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
3005 fpt.SetParameters(7.24,0.120);
3007 for(Int_t nt=0; nt<nTracks; nt++){
3008 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
3009 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
3010 Double_t pt = fpt.GetRandom(); // momentum for f1
3011 // printf("phi %lf eta %lf pt %lf\n",phi,eta,pt);
3013 if(gRandom->Rndm() < 0.5){
3019 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
3021 pxyz[0]=pt*TMath::Cos(phi);
3022 pxyz[1]=pt*TMath::Sin(phi);
3023 pxyz[2]=pt*TMath::Tan(theta);
3024 Double_t cv[21]={0};
3025 AliExternalTrackParam *t= new AliExternalTrackParam(orgVertex, pxyz, cv, sign);
3029 AliExternalTrackParam *td = FitDistortedTrack(*t, refX, dir, NULL);
3031 if (pcstream) (*pcstream)<<"track"<<
3037 if(( eta>0.07 )&&( eta<etaCuts )) { // - log(tan(0.5*theta)), theta = 0.5*pi - ATan(5.0/80.0)
3041 Int_t nn=aTrk.GetEntriesFast();
3044 }else if(( eta<-0.07 )&&( eta>-etaCuts )){
3048 Int_t nn=cTrk.GetEntriesFast();
3053 }// end of track loop
3055 vertexer->SetTPCMode();
3056 vertexer->SetConstraintOff();
3058 aV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&daTrk,aId));
3059 avOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&aTrk,aId));
3060 cV = *((AliESDVertex*)vertexer->FindPrimaryVertex(&dcTrk,cId));
3061 cvOrg = *((AliESDVertex*)vertexer->FindPrimaryVertex(&cTrk,cId));
3062 if (pcstream) (*pcstream)<<"vertex"<<
3063 "x="<<orgVertex[0]<<
3064 "y="<<orgVertex[1]<<
3065 "z="<<orgVertex[2]<<
3066 "av.="<<&aV<< // distorted vertex A side
3067 "cv.="<<&cV<< // distroted vertex C side
3068 "avO.="<<&avOrg<< // original vertex A side
3075 void AliTPCCorrection::AddVisualCorrection(AliTPCCorrection* corr, Int_t position){
3077 // make correction available for visualization using
3078 // TFormula, TFX and TTree::Draw
3079 // important in order to check corrections and also compute dervied variables
3080 // e.g correction partial derivatives
3082 // NOTE - class is not owner of correction
3084 if (!fgVisualCorrection) fgVisualCorrection=new TObjArray(10000);
3085 if (position>=fgVisualCorrection->GetEntriesFast())
3086 fgVisualCorrection->Expand((position+10)*2);
3087 fgVisualCorrection->AddAt(corr, position);
3090 AliTPCCorrection* AliTPCCorrection::GetVisualCorrection(Int_t position) {
3092 // Get visula correction registered at index=position
3094 return fgVisualCorrection? (AliTPCCorrection*)fgVisualCorrection->At(position):0;
3099 Double_t AliTPCCorrection::GetCorrSector(Double_t sector, Double_t r, Double_t kZ, Int_t axisType, Int_t corrType){
3101 // calculate the correction at given position - check the geffCorr
3103 // corrType return values
3109 if (!fgVisualCorrection) return 0;
3110 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3111 if (!corr) return 0;
3113 Double_t phi=sector*TMath::Pi()/9.;
3114 Double_t gx = r*TMath::Cos(phi);
3115 Double_t gy = r*TMath::Sin(phi);
3117 Int_t nsector=(gz>=0) ? 0:18;
3121 Float_t distPoint[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
3122 corr->DistortPoint(distPoint, nsector);
3123 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3124 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3125 Double_t phi0=TMath::ATan2(gy,gx);
3126 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3127 if (axisType==0) return r1-r0;
3128 if (axisType==1) return (phi1-phi0)*r0;
3129 if (axisType==2) return distPoint[2]-gz;
3130 if (axisType==3) return (TMath::Cos(phi)*(distPoint[0]-gx)+ TMath::Cos(phi)*(distPoint[1]-gy));
3134 Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
3136 // return correction at given x,y,z
3138 if (!fgVisualCorrection) return 0;
3139 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3140 if (!corr) return 0;
3141 Double_t phi0= TMath::ATan2(gy,gx);
3142 Int_t nsector=(gz>=0) ? 0:18;
3143 Float_t distPoint[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
3144 corr->CorrectPoint(distPoint, nsector);
3145 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3146 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3147 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3148 if (axisType==0) return r1-r0;
3149 if (axisType==1) return (phi1-phi0)*r0;
3150 if (axisType==2) return distPoint[2]-gz;
3154 Double_t AliTPCCorrection::GetCorrXYZDz(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType,Double_t delta){
3156 // return correction at given x,y,z
3158 if (!fgVisualCorrection) return 0;
3159 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3160 if (!corr) return 0;
3161 Double_t phi0= TMath::ATan2(gy,gx);
3162 Int_t nsector=(gz>=0) ? 0:18;
3163 Float_t distPoint[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
3164 Float_t dxyz[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
3166 corr->GetCorrectionDz(distPoint, nsector,dxyz,delta);
3167 distPoint[0]+=dxyz[0];
3168 distPoint[1]+=dxyz[1];
3169 distPoint[2]+=dxyz[2];
3170 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3171 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3172 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3173 if (axisType==0) return r1-r0;
3174 if (axisType==1) return (phi1-phi0)*r0;
3175 if (axisType==2) return distPoint[2]-gz;
3179 Double_t AliTPCCorrection::GetCorrXYZIntegrateZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType,Double_t delta){
3181 // return correction at given x,y,z
3183 if (!fgVisualCorrection) return 0;
3184 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3185 if (!corr) return 0;
3186 Double_t phi0= TMath::ATan2(gy,gx);
3187 Int_t nsector=(gz>=0) ? 0:18;
3188 Float_t distPoint[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
3189 Float_t dxyz[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
3191 corr->GetCorrectionIntegralDz(distPoint, nsector,dxyz,delta);
3192 distPoint[0]+=dxyz[0];
3193 distPoint[1]+=dxyz[1];
3194 distPoint[2]+=dxyz[2];
3195 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3196 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3197 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3198 if (axisType==0) return r1-r0;
3199 if (axisType==1) return (phi1-phi0)*r0;
3200 if (axisType==2) return distPoint[2]-gz;
3205 Double_t AliTPCCorrection::GetDistXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType){
3207 // return correction at given x,y,z
3209 if (!fgVisualCorrection) return 0;
3210 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3211 if (!corr) return 0;
3212 Double_t phi0= TMath::ATan2(gy,gx);
3213 Int_t nsector=(gz>=0) ? 0:18;
3214 Float_t distPoint[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
3215 corr->DistortPoint(distPoint, nsector);
3216 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3217 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3218 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3219 if (axisType==0) return r1-r0;
3220 if (axisType==1) return (phi1-phi0)*r0;
3221 if (axisType==2) return distPoint[2]-gz;
3225 Double_t AliTPCCorrection::GetDistXYZDz(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType,Double_t delta){
3227 // return correction at given x,y,z
3229 if (!fgVisualCorrection) return 0;
3230 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3231 if (!corr) return 0;
3232 Double_t phi0= TMath::ATan2(gy,gx);
3233 Int_t nsector=(gz>=0) ? 0:18;
3234 Float_t distPoint[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
3235 Float_t dxyz[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
3237 corr->GetDistortionDz(distPoint, nsector,dxyz,delta);
3238 distPoint[0]+=dxyz[0];
3239 distPoint[1]+=dxyz[1];
3240 distPoint[2]+=dxyz[2];
3241 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3242 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3243 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3244 if (axisType==0) return r1-r0;
3245 if (axisType==1) return (phi1-phi0)*r0;
3246 if (axisType==2) return distPoint[2]-gz;
3250 Double_t AliTPCCorrection::GetDistXYZIntegrateZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType, Int_t corrType,Double_t delta){
3252 // return correction at given x,y,z
3254 if (!fgVisualCorrection) return 0;
3255 AliTPCCorrection *corr = (AliTPCCorrection*)fgVisualCorrection->At(corrType);
3256 if (!corr) return 0;
3257 Double_t phi0= TMath::ATan2(gy,gx);
3258 Int_t nsector=(gz>=0) ? 0:18;
3259 Float_t distPoint[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
3260 Float_t dxyz[3]={static_cast<Float_t>(gx),static_cast<Float_t>(gy),static_cast<Float_t>(gz)};
3262 corr->GetDistortionIntegralDz(distPoint, nsector,dxyz,delta);
3263 distPoint[0]+=dxyz[0];
3264 distPoint[1]+=dxyz[1];
3265 distPoint[2]+=dxyz[2];
3266 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
3267 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
3268 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
3269 if (axisType==0) return r1-r0;
3270 if (axisType==1) return (phi1-phi0)*r0;
3271 if (axisType==2) return distPoint[2]-gz;
3277 void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){
3279 // Make a laser fit tree for global minimization
3281 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
3282 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
3283 if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
3284 correction->AddVisualCorrection(correction,0); //register correction
3286 // AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
3287 //AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
3289 const Double_t cutErrY=0.05;
3290 const Double_t kSigmaCut=4;
3291 // const Double_t cutErrZ=0.03;
3292 const Double_t kEpsilon=0.00000001;
3293 // const Double_t kMaxDist=1.; // max distance - space correction
3298 AliTPCLaserTrack *ltr=0;
3299 AliTPCLaserTrack::LoadTracks();
3300 tree->SetBranchAddress("dY.",&vecdY);
3301 tree->SetBranchAddress("dZ.",&vecdZ);
3302 tree->SetBranchAddress("eY.",&veceY);
3303 tree->SetBranchAddress("eZ.",&veceZ);
3304 tree->SetBranchAddress("LTr.",<r);
3305 Int_t entries= tree->GetEntries();
3306 TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
3307 Double_t bz=AliTrackerBase::GetBz();
3309 // Double_t globalXYZ[3];
3310 //Double_t globalXYZCorr[3];
3311 for (Int_t ientry=0; ientry<entries; ientry++){
3312 tree->GetEntry(ientry);
3313 if (!ltr->GetVecGX()){
3314 ltr->UpdatePoints();
3319 printf("Entry\t%d\n",ientry);
3320 for (Int_t irow0=0; irow0<158; irow0+=1){
3322 TLinearFitter fitter10(4,"hyp3");
3323 TLinearFitter fitter5(2,"hyp1");
3324 Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0];
3325 if (sector<0) continue;
3326 //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue;
3328 Double_t refX= (*ltr->GetVecLX())[irow0];
3329 Int_t firstRow1 = TMath::Max(irow0-10,0);
3330 Int_t lastRow1 = TMath::Min(irow0+10,158);
3331 Double_t padWidth=(irow0<64)?0.4:0.6;
3332 // make long range fit
3333 for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){
3334 if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
3335 if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
3336 if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
3337 Double_t idealX= (*ltr->GetVecLX())[irow1];
3338 Double_t idealY= (*ltr->GetVecLY())[irow1];
3339 // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
3340 Double_t gx= (*ltr->GetVecGX())[irow1];
3341 Double_t gy= (*ltr->GetVecGY())[irow1];
3342 Double_t gz= (*ltr->GetVecGZ())[irow1];
3343 Double_t measY=(*vecdY)[irow1]+idealY;
3344 Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
3345 // deltaR = R distorted -R ideal
3346 Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
3347 fitter10.AddPoint(xxx,measY,1);
3350 Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
3351 Double_t mean10 =0;// fitter10.GetParameter(0);
3352 Double_t slope10 =0;// fitter10.GetParameter(0);
3353 Double_t cosPart10 = 0;// fitter10.GetParameter(2);
3354 Double_t sinPart10 =0;// fitter10.GetParameter(3);
3356 if (fitter10.GetNpoints()>10){
3358 rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
3359 mean10 = fitter10.GetParameter(0);
3360 slope10 = fitter10.GetParameter(1);
3361 cosPart10 = fitter10.GetParameter(2);
3362 sinPart10 = fitter10.GetParameter(3);
3364 // make short range fit
3366 for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){
3367 if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
3368 if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
3369 if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
3370 Double_t idealX= (*ltr->GetVecLX())[irow1];
3371 Double_t idealY= (*ltr->GetVecLY())[irow1];
3372 // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
3373 Double_t gx= (*ltr->GetVecGX())[irow1];
3374 Double_t gy= (*ltr->GetVecGY())[irow1];
3375 Double_t gz= (*ltr->GetVecGZ())[irow1];
3376 Double_t measY=(*vecdY)[irow1]+idealY;
3377 Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
3378 // deltaR = R distorted -R ideal
3379 Double_t expY= mean10+slope10*(idealX+deltaR-refX);
3380 if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue;
3382 Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
3383 Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
3384 fitter5.AddPoint(xxx,measY-corr,1);
3389 if (fitter5.GetNpoints()<8) isOK=kFALSE;
3391 Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
3392 Double_t offset5 =0;// fitter5.GetParameter(0);
3393 Double_t slope5 =0;// fitter5.GetParameter(0);
3396 rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
3397 offset5 = fitter5.GetParameter(0);
3398 slope5 = fitter5.GetParameter(0);
3403 Double_t phi =(*ltr->GetVecPhi())[irow0];
3404 Double_t theta =ltr->GetTgl();
3405 Double_t mean=(vecdY)->GetMatrixArray()[irow0];
3406 Double_t gx=0,gy=0,gz=0;
3407 Double_t snp = (*ltr->GetVecP2())[irow0];
3408 Int_t bundle= ltr->GetBundle();
3409 Int_t id= ltr->GetId();
3410 // Double_t rms = err->GetMatrixArray()[irow];
3412 gx = (*ltr->GetVecGX())[irow0];
3413 gy = (*ltr->GetVecGY())[irow0];
3414 gz = (*ltr->GetVecGZ())[irow0];
3415 Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0);
3416 fitter10.GetParameters(fit10);
3417 fitter5.GetParameters(fit5);
3418 Double_t idealY= (*ltr->GetVecLY())[irow0];
3419 Double_t measY=(*vecdY)[irow0]+idealY;
3420 Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
3421 if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE;
3423 (*pcstream)<<"fitFull"<< // dumpe also intermediate results
3424 "bz="<<bz<< // magnetic filed used
3425 "dtype="<<dtype<< // detector match type
3426 "ptype="<<ptype<< // parameter type
3427 "theta="<<theta<< // theta
3428 "phi="<<phi<< // phi
3429 "snp="<<snp<< // snp
3432 // // "dsec="<<dsec<<
3433 "refX="<<refX<< // reference radius
3434 "gx="<<gx<< // global position
3435 "gy="<<gy<< // global position
3436 "gz="<<gz<< // global position
3437 "dRrec="<<dRrec<< // delta Radius in reconstruction
3438 "id="<<id<< //bundle