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