X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FTPCbase%2FAliTPCCorrection.cxx;fp=TPC%2FTPCbase%2FAliTPCCorrection.cxx;h=fc7d1200ac84587bca7c43a2792e8f3c888acfa3;hb=73ba6874383113e781fdb75dcac129389f052e4e;hp=0000000000000000000000000000000000000000;hpb=15b503ebe89c2447cc0123de3b265ddb9c871e67;p=u%2Fmrichter%2FAliRoot.git
diff --git a/TPC/TPCbase/AliTPCCorrection.cxx b/TPC/TPCbase/AliTPCCorrection.cxx
new file mode 100644
index 00000000000..fc7d1200ac8
--- /dev/null
+++ b/TPC/TPCbase/AliTPCCorrection.cxx
@@ -0,0 +1,3452 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+// _________________________________________________________________
+//
+// Begin_Html
+//
AliTPCCorrection class
+//
+// The AliTPCCorrection class provides a general framework to deal with space point distortions.
+// An correction class which inherits from here is for example AliTPCExBBShape or AliTPCExBTwist.
+// General virtual functions are (for example) CorrectPoint(x,roc) where x is the vector of initial
+// positions in cartesian coordinates and roc represents the read-out chamber number according to
+// the offline numbering convention. The vector x is overwritten with the corrected coordinates.
+// An alternative usage would be CorrectPoint(x,roc,dx), which leaves the vector x untouched, but
+// returns the distortions via the vector dx.
+// This class is normally used via the general class AliTPCComposedCorrection.
+//
+// Furthermore, the class contains basic geometrical descriptions like field cage radii
+// (fgkIFCRadius, fgkOFCRadius) and length (fgkTPCZ0) plus the voltages. Also, the definitions
+// of size and widths of the fulcrums building the grid of the final look-up table, which is
+// then interpolated, is defined in kNX and fgkXList).
+//
+// All physics-model classes below are derived from this class in order to not duplicate code
+// and to allow a uniform treatment of all physics models.
+//
+//
Poisson solver
+// A numerical solver of the Poisson equation (relaxation technique) is implemented for 2-dimensional
+// geometries (r,z) as well as for 3-dimensional problems (r,$\phi$,z). The corresponding function
+// names are PoissonRelaxation?D. The relevant function arguments are the arrays of the boundary and
+// initial conditions (ArrayofArrayV, ArrayofChargeDensities) as well as the grid granularity which
+// is used during the calculation. These inputs can be chosen according to the needs of the physical
+// effect which is supposed to be simulated. In the 3D version, different symmetry conditions can be set
+// in order to reduce the calculation time (used in AliTPCFCVoltError3D).
+//
+//
Unified plotting functionality
+// Generic plot functions were implemented. They return a histogram pointer in the chosen plane of
+// the TPC drift volume with a selectable grid granularity and the magnitude of the correction vector.
+// For example, the function CreateHistoDZinXY(z,nx,ny) returns a 2-dimensional histogram which contains
+// the longitudinal corrections $dz$ in the (x,y)-plane at the given z position with the granularity of
+// nx and ny. The magnitude of the corrections is defined by the class from which this function is called.
+// In the same manner, standard plots for the (r,$\phi$)-plane and for the other corrections like $dr$ and $rd\phi$ are available
+//
+// Note: This class is normally used via the class AliTPCComposedCorrection
+// End_Html
+//
+// Begin_Macro(source)
+// {
+// gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
+// TCanvas *c2 = new TCanvas("cAliTPCCorrection","cAliTPCCorrection",700,1050); c2->Divide(2,3);
+// AliTPCROCVoltError3D roc; // EXAMPLE PLOTS - SEE BELOW
+// roc.SetROCDataFileName("$ALICE_ROOT/TPC/Calib/maps/TPCROCdzSurvey.root");
+// roc.SetOmegaTauT1T2(0,1,1); // B=0
+// Float_t z0 = 1; // at +1 cm -> A side
+// c2->cd(1); roc.CreateHistoDRinXY(1.,300,300)->Draw("cont4z");
+// c2->cd(3);roc.CreateHistoDRPhiinXY(1.,300,300)->Draw("cont4z");
+// c2->cd(5);roc.CreateHistoDZinXY(1.,300,300)->Draw("cont4z");
+// Float_t phi0=0.5;
+// c2->cd(2);roc.CreateHistoDRinZR(phi0)->Draw("surf2");
+// c2->cd(4);roc.CreateHistoDRPhiinZR(phi0)->Draw("surf2");
+// c2->cd(6);roc.CreateHistoDZinZR(phi0)->Draw("surf2");
+// return c2;
+// }
+// End_Macro
+//
+// Begin_Html
+//
+// Date: 27/04/2010
+// Authors: Magnus Mager, Stefan Rossegger, Jim Thomas
+// End_Html
+// _________________________________________________________________
+
+
+#include "Riostream.h"
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include "TVectorD.h"
+#include "AliTPCParamSR.h"
+
+#include "AliTPCCorrection.h"
+#include "AliLog.h"
+
+#include "AliExternalTrackParam.h"
+#include "AliTrackPointArray.h"
+#include "TDatabasePDG.h"
+#include "AliTrackerBase.h"
+#include "AliTPCROC.h"
+#include "THnSparse.h"
+
+#include "AliTPCLaserTrack.h"
+#include "AliESDVertex.h"
+#include "AliVertexerTracks.h"
+#include "TDatabasePDG.h"
+#include "TF1.h"
+#include "TRandom.h"
+
+#include "TDatabasePDG.h"
+
+#include "AliTPCTransform.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCExB.h"
+
+#include "AliTPCRecoParam.h"
+#include "TLinearFitter.h"
+#include
+
+ClassImp(AliTPCCorrection)
+
+
+TObjArray *AliTPCCorrection::fgVisualCorrection=0;
+// instance of correction for visualization
+
+
+// FIXME: the following values should come from the database
+const Double_t AliTPCCorrection::fgkTPCZ0 = 249.7; // nominal gating grid position
+const Double_t AliTPCCorrection::fgkIFCRadius= 83.5; // radius which renders the "18 rod manifold" best -> compare calc. of Jim Thomas
+// compare gkIFCRadius= 83.05: Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
+const Double_t AliTPCCorrection::fgkOFCRadius= 254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
+const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
+const Double_t AliTPCCorrection::fgkCathodeV = -100000.0; // Cathode Voltage (volts)
+const Double_t AliTPCCorrection::fgkGG = -70.0; // Gating Grid voltage (volts)
+
+const Double_t AliTPCCorrection::fgkdvdE = 0.0024; // [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)
+
+const Double_t AliTPCCorrection::fgkEM = -1.602176487e-19/9.10938215e-31; // charge/mass in [C/kg]
+const Double_t AliTPCCorrection::fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)]
+
+
+AliTPCCorrection::AliTPCCorrection()
+ : TNamed("correction_unity","unity"),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1), fIsLocal(kFALSE)
+{
+ //
+ // default constructor
+ //
+ if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
+
+ InitLookUpfulcrums();
+
+}
+
+AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
+ : TNamed(name,title),fILow(0),fJLow(0),fKLow(0), fT1(1), fT2(1), fIsLocal(kFALSE)
+{
+ //
+ // default constructor, that set the name and title
+ //
+ if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
+
+ InitLookUpfulcrums();
+
+}
+
+AliTPCCorrection::~AliTPCCorrection() {
+ //
+ // virtual destructor
+ //
+}
+
+Bool_t AliTPCCorrection::AddCorrectionCompact(AliTPCCorrection* corr, Double_t weight){
+ //
+ // Add correction and make them compact
+ // Assumptions:
+ // - origin of distortion/correction are additive
+ // - only correction ot the same type supported ()
+ if (corr==NULL) {
+ AliError("Zerro pointer - correction");
+ return kFALSE;
+ }
+ AliError(TString::Format("Correction %s not implementend",IsA()->GetName()).Data());
+ return kFALSE;
+}
+
+
+void AliTPCCorrection::CorrectPoint(Float_t x[], Short_t roc) {
+ //
+ // Corrects the initial coordinates x (cartesian coordinates)
+ // according to the given effect (inherited classes)
+ // roc represents the TPC read out chamber (offline numbering convention)
+ //
+ Float_t dx[3];
+ GetCorrection(x,roc,dx);
+ for (Int_t j=0;j<3;++j) x[j]+=dx[j];
+}
+
+void AliTPCCorrection::CorrectPoint(const Float_t x[], Short_t roc,Float_t xp[]) {
+ //
+ // Corrects the initial coordinates x (cartesian coordinates) and stores the new
+ // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
+ // roc represents the TPC read out chamber (offline numbering convention)
+ //
+ Float_t dx[3];
+ GetCorrection(x,roc,dx);
+ for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
+}
+
+void AliTPCCorrection::DistortPoint(Float_t x[], Short_t roc) {
+ //
+ // Distorts the initial coordinates x (cartesian coordinates)
+ // according to the given effect (inherited classes)
+ // roc represents the TPC read out chamber (offline numbering convention)
+ //
+ Float_t dx[3];
+ GetDistortion(x,roc,dx);
+ for (Int_t j=0;j<3;++j) x[j]+=dx[j];
+}
+
+void AliTPCCorrection::DistortPointLocal(Float_t x[], Short_t roc) {
+ //
+ // Distorts the initial coordinates x (cartesian coordinates)
+ // according to the given effect (inherited classes)
+ // roc represents the TPC read out chamber (offline numbering convention)
+ //
+ Float_t gxyz[3]={0,0,0};
+ Double_t alpha = TMath::TwoPi()*(roc%18+0.5)/18;
+ Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
+ gxyz[0]= ca*x[0]+sa*x[1];
+ gxyz[1]= -sa*x[0]+ca*x[1];
+ gxyz[2]= x[2];
+ DistortPoint(gxyz,roc);
+ x[0]= ca*gxyz[0]-sa*gxyz[1];
+ x[1]= +sa*gxyz[0]+ca*gxyz[1];
+ x[2]= gxyz[2];
+}
+void AliTPCCorrection::CorrectPointLocal(Float_t x[], Short_t roc) {
+ //
+ // Distorts the initial coordinates x (cartesian coordinates)
+ // according to the given effect (inherited classes)
+ // roc represents the TPC read out chamber (offline numbering convention)
+ //
+ Float_t gxyz[3]={0,0,0};
+ Double_t alpha = TMath::TwoPi()*(roc%18+0.5)/18;
+ Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
+ gxyz[0]= ca*x[0]+sa*x[1];
+ gxyz[1]= -sa*x[0]+ca*x[1];
+ gxyz[2]= x[2];
+ CorrectPoint(gxyz,roc);
+ x[0]= ca*gxyz[0]-sa*gxyz[1];
+ x[1]= sa*gxyz[0]+ca*gxyz[1];
+ x[2]= gxyz[2];
+}
+
+void AliTPCCorrection::DistortPoint(const Float_t x[], Short_t roc,Float_t xp[]) {
+ //
+ // Distorts the initial coordinates x (cartesian coordinates) and stores the new
+ // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
+ // roc represents the TPC read out chamber (offline numbering convention)
+ //
+ Float_t dx[3];
+ GetDistortion(x,roc,dx);
+ for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
+}
+
+void AliTPCCorrection::GetCorrection(const Float_t /*x*/[], Short_t /*roc*/,Float_t dx[]) {
+ //
+ // This function delivers the correction values dx in respect to the inital coordinates x
+ // roc represents the TPC read out chamber (offline numbering convention)
+ // Note: The dx is overwritten by the inherited effectice class ...
+ //
+ for (Int_t j=0;j<3;++j) { dx[j]=0.; }
+}
+
+void AliTPCCorrection::GetDistortion(const Float_t x[], Short_t roc,Float_t dx[]) {
+ //
+ // This function delivers the distortion values dx in respect to the inital coordinates x
+ // roc represents the TPC read out chamber (offline numbering convention)
+ //
+ GetCorrection(x,roc,dx);
+ for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
+}
+
+void AliTPCCorrection::GetCorrectionDz(const Float_t x[], Short_t roc,Float_t dx[], Float_t delta) {
+ // author: marian.ivanov@cern.ch
+ //
+ // In this (virtual)function calculates the dx'/dz, dy'/dz and dz'/dz at given point (x,y,z)
+ // Generic implementation. Better precision can be acchieved knowing the internal structure
+ // of underlying trasnformation. Derived classes can reimplement it.
+ // To calculate correction is fitted in small neighberhood:
+ // (x+-delta,y+-delta,z+-delta) where delta is an argument
+ //
+ // Input parameters:
+ // x[] - space point corrdinate
+ // roc - readout chamber identifier (important e.g to do not miss the side of detector)
+ // delta - define the size of neighberhood
+ // Output parameter:
+ // dx[] - array {dx'/dz, dy'/dz , dz'/dz }
+
+ // if (fIsLocal){ //standard implemenation provides the correction/distortion integrated over full drift length
+ //
+ //
+ // GetCorrection(xyz,roc,dxyz);
+ // }
+ static TLinearFitter fitx(2,"pol1");
+ static TLinearFitter fity(2,"pol1");
+ static TLinearFitter fitz(2,"pol1");
+ fitx.ClearPoints();
+ fity.ClearPoints();
+ fitz.ClearPoints();
+ Int_t zmin=-2;
+ Int_t zmax=0;
+ //adjust limits around CE to stay on one side
+ if ((roc%36)<18) {
+ //A-Side
+ if ((x[2]+zmin*delta)<0){
+ zmin=0;
+ zmax=2;
+ if ((x[2]-delta)>0){
+ zmin=-1;
+ zmax=1;
+ }
+ }
+ } else {
+ //C-Side
+ zmin=0;
+ zmax=2;
+ if ((x[2]+zmax*delta)>0){
+ zmin=-2;
+ zmax=0;
+ if ((x[2]+delta)<0){
+ zmin=-1;
+ zmax=1;
+ }
+ }
+ }
+
+ for (Int_t xdelta=-1; xdelta<=1; xdelta++)
+ for (Int_t ydelta=-1; ydelta<=1; ydelta++){
+// for (Int_t zdelta=-1; zdelta<=1; zdelta++){
+// for (Int_t xdelta=-2; xdelta<=0; xdelta++)
+// for (Int_t ydelta=-2; ydelta<=0; ydelta++){
+ for (Int_t zdelta=zmin; zdelta<=zmax; zdelta++){
+ //TODO: what happens if x[2] is on the A-Side, but x[2]+zdelta*delta
+ // will be on the C-Side?
+ Float_t xyz[3]={x[0]+xdelta*delta, x[1]+ydelta*delta, x[2]+zdelta*delta};
+ Float_t dxyz[3];
+ GetCorrection(xyz,roc,dxyz);
+ Double_t adelta=zdelta*delta;
+ fitx.AddPoint(&adelta, dxyz[0]);
+ fity.AddPoint(&adelta, dxyz[1]);
+ fitz.AddPoint(&adelta, dxyz[2]);
+ }
+ }
+ fitx.Eval();
+ fity.Eval();
+ fitz.Eval();
+ dx[0] = fitx.GetParameter(1);
+ dx[1] = fity.GetParameter(1);
+ dx[2] = fitz.GetParameter(1);
+}
+
+void AliTPCCorrection::GetDistortionDz(const Float_t x[], Short_t roc,Float_t dx[], Float_t delta) {
+ // author: marian.ivanov@cern.ch
+ //
+ // In this (virtual)function calculates the dx'/dz, dy'/dz and dz'/dz at given point (x,y,z)
+ // Generic implementation. Better precision can be acchieved knowing the internal structure
+ // of underlying trasnformation. Derived classes can reimplement it.
+ // To calculate distortion is fitted in small neighberhood:
+ // (x+-delta,y+-delta,z+-delta) where delta is an argument
+ //
+ // Input parameters:
+ // x[] - space point corrdinate
+ // roc - readout chamber identifier (important e.g to do not miss the side of detector)
+ // delta - define the size of neighberhood
+ // Output parameter:
+ // dx[] - array {dx'/dz, dy'/dz , dz'/dz }
+
+ static TLinearFitter fitx(2,"pol1");
+ static TLinearFitter fity(2,"pol1");
+ static TLinearFitter fitz(2,"pol1");
+ fitx.ClearPoints();
+ fity.ClearPoints();
+ fitz.ClearPoints();
+
+ Int_t zmin=-1;
+ Int_t zmax=1;
+ //adjust limits around CE to stay on one side
+ if ((roc%36)<18) {
+ //A-Side
+ if ((x[2]+zmin*delta)<0){
+ zmin=0;
+ zmax=2;
+ }
+ } else {
+ //C-Side
+ if ((x[2]+zmax*delta)>0){
+ zmin=-2;
+ zmax=0;
+ }
+ }
+
+ //TODO: in principle one shuld check that x[2]+zdelta*delta does not get 'out of' bounds,
+ // so close to the CE it doesn't change the sign, since then the corrections will be wrong ...
+ for (Int_t xdelta=-1; xdelta<=1; xdelta++)
+ for (Int_t ydelta=-1; ydelta<=1; ydelta++){
+ for (Int_t zdelta=zmin; zdelta<=zmax; zdelta++){
+ //TODO: what happens if x[2] is on the A-Side, but x[2]+zdelta*delta
+ // will be on the C-Side?
+ //TODO: For the C-Side, does this have the correct sign?
+ Float_t xyz[3]={x[0]+xdelta*delta, x[1]+ydelta*delta, x[2]+zdelta*delta};
+ Float_t dxyz[3];
+ GetDistortion(xyz,roc,dxyz);
+ Double_t adelta=zdelta*delta;
+ fitx.AddPoint(&adelta, dxyz[0]);
+ fity.AddPoint(&adelta, dxyz[1]);
+ fitz.AddPoint(&adelta, dxyz[2]);
+ }
+ }
+ fitx.Eval();
+ fity.Eval();
+ fitz.Eval();
+ dx[0] = fitx.GetParameter(1);
+ dx[1] = fity.GetParameter(1);
+ dx[2] = fitz.GetParameter(1);
+}
+
+void AliTPCCorrection::GetCorrectionIntegralDz(const Float_t x[], Short_t roc,Float_t dx[], Float_t delta){
+ //
+ // Integrate 3D distortion along drift lines starting from the roc plane
+ // to the expected z position of the point, this assumes that dz is small
+ // and the error propagating to z' instead of the correct z is negligible
+ // To define the drift lines virtual function AliTPCCorrection::GetCorrectionDz is used
+ //
+ // Input parameters:
+ // x[] - space point corrdinate
+ // roc - readout chamber identifier (important e.g to do not miss the side of detector)
+ // delta - define the size of neighberhood
+ // Output parameter:
+ // dx[] - array { integral(dx'/dz), integral(dy'/dz) , integral(dz'/dz) }
+
+ Float_t zroc= ((roc%36)<18) ? fgkTPCZ0:-fgkTPCZ0;
+ Double_t zdrift = TMath::Abs(x[2]-zroc);
+ Int_t nsteps = Int_t(zdrift/delta)+1;
+ //
+ //
+ Float_t xyz[3]={x[0],x[1],zroc};
+ Float_t dxyz[3]={x[0],x[1],x[2]};
+ Short_t side=(roc/18)%2;
+ Float_t sign=1-2*side;
+ Double_t sumdz=0;
+ for (Int_t i=0;ifgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-fgkTPCZ0);
+// if (xyz[2]-deltaZ<-fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-fgkTPCZ0);
+ // protect again integrating through the CE
+ if (side==0){
+ if (xyz[2]+deltaZ<0) deltaZ=-xyz[2]+1e-20;
+ } else {
+ if (xyz[2]+deltaZ>0) deltaZ=xyz[2]-+1e-20;
+ }
+ // since at larger drift (smaller z) the corrections are larger (absolute, but negative)
+ // the slopes will be positive.
+ // but since we chose deltaZ opposite sign the singn of the corretion should be fine
+
+ Float_t xyz2[3]={xyz[0],xyz[1],static_cast(xyz[2]+deltaZ/2.)};
+ GetCorrectionDz(xyz2,roc,dxyz,delta/2.);
+ xyz[0]+=deltaZ*dxyz[0];
+ xyz[1]+=deltaZ*dxyz[1];
+ xyz[2]+=deltaZ; //
+ sumdz+=deltaZ*dxyz[2];
+ }
+ //
+ dx[0]=xyz[0]-x[0];
+ dx[1]=xyz[1]-x[1];
+ dx[2]= sumdz; //TODO: is sumdz correct?
+}
+
+void AliTPCCorrection::GetDistortionIntegralDz(const Float_t x[], Short_t roc,Float_t dx[], Float_t delta){
+ //
+ // Integrate 3D distortion along drift lines
+ // To define the drift lines virtual function AliTPCCorrection::GetCorrectionDz is used
+ //
+ // Input parameters:
+ // x[] - space point corrdinate
+ // roc - readout chamber identifier (important e.g to do not miss the side of detector)
+ // delta - define the size of neighberhood
+ // Output parameter:
+ // dx[] - array { integral(dx'/dz), integral(dy'/dz) , integral(dz'/dz) }
+
+ Float_t zroc= ((roc%36)<18) ? fgkTPCZ0:-fgkTPCZ0;
+ Double_t zdrift = TMath::Abs(x[2]-zroc);
+ Int_t nsteps = Int_t(zdrift/delta)+1;
+ //
+ //
+ Float_t xyz[3]={x[0],x[1],x[2]};
+ Float_t dxyz[3]={x[0],x[1],x[2]};
+ Float_t sign=((roc%36)<18) ? 1.:-1.;
+ Double_t sumdz=0;
+ for (Int_t i=0;ifgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-zroc);
+ if (xyz[2]-deltaZ<-fgkTPCZ0) deltaZ=TMath::Abs(xyz[2]-zroc);
+ // since at larger drift (smaller z) the distortions are larger
+ // the slopes will be negative.
+ // and since we are moving towards the read-out plane the deltaZ for
+ // weighting the dK/dz should have the opposite sign
+ deltaZ*=sign;
+ Float_t xyz2[3]={xyz[0],xyz[1],static_cast(xyz[2]+deltaZ/2.)};
+ GetDistortionDz(xyz2,roc,dxyz,delta/2.);
+ xyz[0]+=-deltaZ*dxyz[0];
+ xyz[1]+=-deltaZ*dxyz[1];
+ xyz[2]+=deltaZ; //TODO: Should this also be corrected for the dxyz[2]
+ sumdz+=-deltaZ*dxyz[2];
+ }
+ //
+ dx[0]=xyz[0]-x[0];
+ dx[1]=xyz[1]-x[1];
+ dx[2]= sumdz; //TODO: is sumdz correct?
+
+}
+
+
+void AliTPCCorrection::Init() {
+ //
+ // Initialization funtion (not used at the moment)
+ //
+}
+
+void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
+ //
+ // Update function
+ //
+}
+
+void AliTPCCorrection::Print(Option_t* /*option*/) const {
+ //
+ // Print function to check which correction classes are used
+ // option=="d" prints details regarding the setted magnitude
+ // option=="a" prints the C0 and C1 coefficents for calibration purposes
+ //
+ printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
+}
+
+void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
+ //
+ // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
+ // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
+ // calibration run
+ //
+ fT1=t1;
+ fT2=t2;
+ //SetOmegaTauT1T2(omegaTau, t1, t2);
+}
+
+TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
+ //
+ // Simple plot functionality.
+ // Returns a 2d hisogram which represents the corrections in radial direction (dr)
+ // in respect to position z within the XY plane.
+ // The histogramm has nx times ny entries.
+ //
+ AliTPCParam* tpcparam = new AliTPCParamSR;
+
+ TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
+ nx,-250.,250.,ny,-250.,250.);
+ Float_t x[3],dx[3];
+ x[2]=z;
+ Int_t roc=z>0.?0:18; // FIXME
+ for (Int_t iy=1;iy<=ny;++iy) {
+ x[1]=h->GetYaxis()->GetBinCenter(iy);
+ for (Int_t ix=1;ix<=nx;++ix) {
+ x[0]=h->GetXaxis()->GetBinCenter(ix);
+ GetCorrection(x,roc,dx);
+ Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
+ if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
+ Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
+ h->SetBinContent(ix,iy,r1-r0);
+ }
+ else
+ h->SetBinContent(ix,iy,0.);
+ }
+ }
+ delete tpcparam;
+ return h;
+}
+
+TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
+ //
+ // Simple plot functionality.
+ // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
+ // in respect to position z within the XY plane.
+ // The histogramm has nx times ny entries.
+ //
+
+ AliTPCParam* tpcparam = new AliTPCParamSR;
+
+ TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
+ nx,-250.,250.,ny,-250.,250.);
+ Float_t x[3],dx[3];
+ x[2]=z;
+ Int_t roc=z>0.?0:18; // FIXME
+ for (Int_t iy=1;iy<=ny;++iy) {
+ x[1]=h->GetYaxis()->GetBinCenter(iy);
+ for (Int_t ix=1;ix<=nx;++ix) {
+ x[0]=h->GetXaxis()->GetBinCenter(ix);
+ GetCorrection(x,roc,dx);
+ Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
+ if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
+ Float_t phi0=TMath::ATan2(x[1] ,x[0] );
+ Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
+
+ Float_t dphi=phi1-phi0;
+ if (dphiTMath::Pi()) dphi-=TMath::TwoPi();
+
+ h->SetBinContent(ix,iy,r0*dphi);
+ }
+ else
+ h->SetBinContent(ix,iy,0.);
+ }
+ }
+ delete tpcparam;
+ return h;
+}
+
+TH2F* AliTPCCorrection::CreateHistoDZinXY(Float_t z,Int_t nx,Int_t ny) {
+ //
+ // Simple plot functionality.
+ // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
+ // in respect to position z within the XY plane.
+ // The histogramm has nx times ny entries.
+ //
+
+ AliTPCParam* tpcparam = new AliTPCParamSR;
+
+ TH2F *h=CreateTH2F("dz_xy",GetTitle(),"x [cm]","y [cm]","dz [cm]",
+ nx,-250.,250.,ny,-250.,250.);
+ Float_t x[3],dx[3];
+ x[2]=z;
+ Int_t roc=z>0.?0:18; // FIXME
+ for (Int_t iy=1;iy<=ny;++iy) {
+ x[1]=h->GetYaxis()->GetBinCenter(iy);
+ for (Int_t ix=1;ix<=nx;++ix) {
+ x[0]=h->GetXaxis()->GetBinCenter(ix);
+ GetCorrection(x,roc,dx);
+ Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
+ if (tpcparam->GetPadRowRadii(0,0)<=r0 && r0<=tpcparam->GetPadRowRadii(36,95)) {
+ h->SetBinContent(ix,iy,dx[2]);
+ }
+ else
+ h->SetBinContent(ix,iy,0.);
+ }
+ }
+ delete tpcparam;
+ return h;
+}
+
+TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
+ //
+ // Simple plot functionality.
+ // Returns a 2d hisogram which represents the corrections in r direction (dr)
+ // in respect to angle phi within the ZR plane.
+ // The histogramm has nx times ny entries.
+ //
+ TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
+ nz,-250.,250.,nr,85.,250.);
+ Float_t x[3],dx[3];
+ for (Int_t ir=1;ir<=nr;++ir) {
+ Float_t radius=h->GetYaxis()->GetBinCenter(ir);
+ x[0]=radius*TMath::Cos(phi);
+ x[1]=radius*TMath::Sin(phi);
+ for (Int_t iz=1;iz<=nz;++iz) {
+ x[2]=h->GetXaxis()->GetBinCenter(iz);
+ Int_t roc=x[2]>0.?0:18; // FIXME
+ GetCorrection(x,roc,dx);
+ Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
+ Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
+ h->SetBinContent(iz,ir,r1-r0);
+ }
+ }
+ return h;
+
+}
+
+TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
+ //
+ // Simple plot functionality.
+ // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
+ // in respect to angle phi within the ZR plane.
+ // The histogramm has nx times ny entries.
+ //
+ TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
+ nz,-250.,250.,nr,85.,250.);
+ Float_t x[3],dx[3];
+ for (Int_t iz=1;iz<=nz;++iz) {
+ x[2]=h->GetXaxis()->GetBinCenter(iz);
+ Int_t roc=x[2]>0.?0:18; // FIXME
+ for (Int_t ir=1;ir<=nr;++ir) {
+ Float_t radius=h->GetYaxis()->GetBinCenter(ir);
+ x[0]=radius*TMath::Cos(phi);
+ x[1]=radius*TMath::Sin(phi);
+ GetCorrection(x,roc,dx);
+ Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
+ Float_t phi0=TMath::ATan2(x[1] ,x[0] );
+ Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
+
+ Float_t dphi=phi1-phi0;
+ if (dphiTMath::Pi()) dphi-=TMath::TwoPi();
+
+ h->SetBinContent(iz,ir,r0*dphi);
+ }
+ }
+ return h;
+}
+
+TH2F* AliTPCCorrection::CreateHistoDZinZR(Float_t phi,Int_t nz,Int_t nr) {
+ //
+ // Simple plot functionality.
+ // Returns a 2d hisogram which represents the corrections in longitudinal direction (dz)
+ // in respect to angle phi within the ZR plane.
+ // The histogramm has nx times ny entries.
+ //
+ TH2F *h=CreateTH2F("dz_zr",GetTitle(),"z [cm]","r [cm]","dz [cm]",
+ nz,-250.,250.,nr,85.,250.);
+ Float_t x[3],dx[3];
+ for (Int_t ir=1;ir<=nr;++ir) {
+ Float_t radius=h->GetYaxis()->GetBinCenter(ir);
+ x[0]=radius*TMath::Cos(phi);
+ x[1]=radius*TMath::Sin(phi);
+ for (Int_t iz=1;iz<=nz;++iz) {
+ x[2]=h->GetXaxis()->GetBinCenter(iz);
+ Int_t roc=x[2]>0.?0:18; // FIXME
+ GetCorrection(x,roc,dx);
+ h->SetBinContent(iz,ir,dx[2]);
+ }
+ }
+ return h;
+
+}
+
+
+TH2F* AliTPCCorrection::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) {
+ //
+ // Helper function to create a 2d histogramm of given size
+ //
+
+ TString hname=name;
+ Int_t i=0;
+ if (gDirectory) {
+ while (gDirectory->FindObject(hname.Data())) {
+ hname =name;
+ hname+="_";
+ hname+=i;
+ ++i;
+ }
+ }
+ TH2F *h=new TH2F(hname.Data(),title,
+ nbinsx,xlow,xup,
+ nbinsy,ylow,yup);
+ h->GetXaxis()->SetTitle(xlabel);
+ h->GetYaxis()->SetTitle(ylabel);
+ h->GetZaxis()->SetTitle(zlabel);
+ h->SetStats(0);
+ return h;
+}
+
+// Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
+
+void AliTPCCorrection::Interpolate2DEdistortion( Int_t order, Double_t r, Double_t z,
+ const Double_t er[kNZ][kNR], Double_t &erValue ) {
+ //
+ // Interpolate table - 2D interpolation
+ //
+ Double_t saveEr[5] = {0,0,0,0,0};
+
+ Search( kNZ, fgkZList, z, fJLow ) ;
+ Search( kNR, fgkRList, r, fKLow ) ;
+ if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
+ if ( fKLow < 0 ) fKLow = 0 ;
+ if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
+ if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
+
+ for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
+ saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
+ }
+ erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
+
+}
+
+void AliTPCCorrection::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) {
+ //
+ // Interpolate table - 3D interpolation
+ //
+
+ Double_t saveEr[5]= {0,0,0,0,0};
+ Double_t savedEr[5]= {0,0,0,0,0} ;
+
+ Double_t saveEphi[5]= {0,0,0,0,0};
+ Double_t savedEphi[5]= {0,0,0,0,0} ;
+
+ Double_t saveEz[5]= {0,0,0,0,0};
+ Double_t savedEz[5]= {0,0,0,0,0} ;
+
+ Search( kNZ, fgkZList, z, fILow ) ;
+ Search( kNPhi, fgkPhiList, z, fJLow ) ;
+ Search( kNR, fgkRList, r, fKLow ) ;
+
+ if ( fILow < 0 ) fILow = 0 ; // check if out of range
+ if ( fJLow < 0 ) fJLow = 0 ;
+ if ( fKLow < 0 ) fKLow = 0 ;
+
+ if ( fILow + order >= kNZ - 1 ) fILow = kNZ - 1 - order ;
+ if ( fJLow + order >= kNPhi - 1 ) fJLow = kNPhi - 1 - order ;
+ if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
+
+ for ( Int_t i = fILow ; i < fILow + order + 1 ; i++ ) {
+ for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
+ saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[i][j][fKLow], order, r ) ;
+ saveEphi[j-fJLow] = Interpolate( &fgkRList[fKLow], &ephi[i][j][fKLow], order, r ) ;
+ saveEz[j-fJLow] = Interpolate( &fgkRList[fKLow], &ez[i][j][fKLow], order, r ) ;
+ }
+ savedEr[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEr, order, phi ) ;
+ savedEphi[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEphi, order, phi ) ;
+ savedEz[i-fILow] = Interpolate( &fgkPhiList[fJLow], saveEz, order, phi ) ;
+ }
+ erValue = Interpolate( &fgkZList[fILow], savedEr, order, z ) ;
+ ephiValue = Interpolate( &fgkZList[fILow], savedEphi, order, z ) ;
+ ezValue = Interpolate( &fgkZList[fILow], savedEz, order, z ) ;
+
+}
+
+Double_t AliTPCCorrection::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 ) {
+ //
+ // Interpolate table (TMatrix format) - 2D interpolation
+ //
+
+ static Int_t jlow = 0, klow = 0 ;
+ Double_t saveArray[5] = {0,0,0,0,0} ;
+
+ Search( nx, xv, x, jlow ) ;
+ Search( ny, yv, y, klow ) ;
+ if ( jlow < 0 ) jlow = 0 ; // check if out of range
+ if ( klow < 0 ) klow = 0 ;
+ if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
+ if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
+
+ for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
+ {
+ Double_t *ajkl = &((TMatrixD&)array)(j,klow);
+ saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
+ }
+
+ return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
+
+}
+
+Double_t AliTPCCorrection::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 ) {
+ //
+ // Interpolate table (TMatrix format) - 3D interpolation
+ //
+
+ static Int_t ilow = 0, jlow = 0, klow = 0 ;
+ Double_t saveArray[5]= {0,0,0,0,0};
+ Double_t savedArray[5]= {0,0,0,0,0} ;
+
+ Search( nx, xv, x, ilow ) ;
+ Search( ny, yv, y, jlow ) ;
+ Search( nz, zv, z, klow ) ;
+
+ if ( ilow < 0 ) ilow = 0 ; // check if out of range
+ if ( jlow < 0 ) jlow = 0 ;
+ if ( klow < 0 ) klow = 0 ;
+
+ if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
+ if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
+ if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
+
+ for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
+ {
+ TMatrixD &table = *arrayofArrays[k] ;
+ for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
+ {
+ saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
+ }
+ savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
+ }
+ return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
+
+}
+
+Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
+ Int_t order, Double_t x ) {
+ //
+ // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
+ //
+
+ Double_t y ;
+ if ( order == 2 ) { // Quadratic Interpolation = 2
+ y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
+ y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
+ y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
+ } else { // Linear Interpolation = 1
+ y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
+ }
+
+ return (y);
+
+}
+
+Float_t AliTPCCorrection::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 TMatrixF &array ) {
+ //
+ // Interpolate table (TMatrix format) - 2D interpolation
+ // Float version (in order to decrease the OCDB size)
+ //
+
+ static Int_t jlow = 0, klow = 0 ;
+ Float_t saveArray[5] = {0.,0.,0.,0.,0.} ;
+
+ Search( nx, xv, x, jlow ) ;
+ Search( ny, yv, y, klow ) ;
+ if ( jlow < 0 ) jlow = 0 ; // check if out of range
+ if ( klow < 0 ) klow = 0 ;
+ if ( jlow + order >= nx - 1 ) jlow = nx - 1 - order ;
+ if ( klow + order >= ny - 1 ) klow = ny - 1 - order ;
+
+ for ( Int_t j = jlow ; j < jlow + order + 1 ; j++ )
+ {
+ Float_t *ajkl = &((TMatrixF&)array)(j,klow);
+ saveArray[j-jlow] = Interpolate( &yv[klow], ajkl , order, y ) ;
+ }
+
+ return( Interpolate( &xv[jlow], saveArray, order, x ) ) ;
+
+}
+
+Float_t AliTPCCorrection::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[],
+ TMatrixF **arrayofArrays ) {
+ //
+ // Interpolate table (TMatrix format) - 3D interpolation
+ // Float version (in order to decrease the OCDB size)
+ //
+
+ static Int_t ilow = 0, jlow = 0, klow = 0 ;
+ Float_t saveArray[5]= {0.,0.,0.,0.,0.};
+ Float_t savedArray[5]= {0.,0.,0.,0.,0.} ;
+
+ Search( nx, xv, x, ilow ) ;
+ Search( ny, yv, y, jlow ) ;
+ Search( nz, zv, z, klow ) ;
+
+ if ( ilow < 0 ) ilow = 0 ; // check if out of range
+ if ( jlow < 0 ) jlow = 0 ;
+ if ( klow < 0 ) klow = 0 ;
+
+ if ( ilow + order >= nx - 1 ) ilow = nx - 1 - order ;
+ if ( jlow + order >= ny - 1 ) jlow = ny - 1 - order ;
+ if ( klow + order >= nz - 1 ) klow = nz - 1 - order ;
+
+ for ( Int_t k = klow ; k < klow + order + 1 ; k++ )
+ {
+ TMatrixF &table = *arrayofArrays[k] ;
+ for ( Int_t i = ilow ; i < ilow + order + 1 ; i++ )
+ {
+ saveArray[i-ilow] = Interpolate( &yv[jlow], &table(i,jlow), order, y ) ;
+ }
+ savedArray[k-klow] = Interpolate( &xv[ilow], saveArray, order, x ) ;
+ }
+ return( Interpolate( &zv[klow], savedArray, order, z ) ) ;
+
+}
+Float_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Float_t yArray[],
+ Int_t order, Double_t x ) {
+ //
+ // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
+ // Float version (in order to decrease the OCDB size)
+ //
+
+ Float_t y ;
+ if ( order == 2 ) { // Quadratic Interpolation = 2
+ y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
+ y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
+ y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
+ } else { // Linear Interpolation = 1
+ y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
+ }
+
+ return (y);
+
+}
+
+
+
+void AliTPCCorrection::Search( Int_t n, const Double_t xArray[], Double_t x, Int_t &low ) {
+ //
+ // Search an ordered table by starting at the most recently used point
+ //
+
+ Long_t middle, high ;
+ Int_t ascend = 0, increment = 1 ;
+
+ if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
+
+ if ( low < 0 || low > n-1 ) {
+ low = -1 ; high = n ;
+ } else { // Ordered Search phase
+ if ( (Int_t)( x >= xArray[low] ) == ascend ) {
+ if ( low == n-1 ) return ;
+ high = low + 1 ;
+ while ( (Int_t)( x >= xArray[high] ) == ascend ) {
+ low = high ;
+ increment *= 2 ;
+ high = low + increment ;
+ if ( high > n-1 ) { high = n ; break ; }
+ }
+ } else {
+ if ( low == 0 ) { low = -1 ; return ; }
+ high = low - 1 ;
+ while ( (Int_t)( x < xArray[low] ) == ascend ) {
+ high = low ;
+ increment *= 2 ;
+ if ( increment >= high ) { low = -1 ; break ; }
+ else low = high - increment ;
+ }
+ }
+ }
+
+ while ( (high-low) != 1 ) { // Binary Search Phase
+ middle = ( high + low ) / 2 ;
+ if ( (Int_t)( x >= xArray[middle] ) == ascend )
+ low = middle ;
+ else
+ high = middle ;
+ }
+
+ if ( x == xArray[n-1] ) low = n-2 ;
+ if ( x == xArray[0] ) low = 0 ;
+
+}
+
+void AliTPCCorrection::InitLookUpfulcrums() {
+ //
+ // Initialization of interpolation points - for main look up table
+ // (course grid in the middle, fine grid on the borders)
+ //
+
+ AliTPCROC * roc = AliTPCROC::Instance();
+ const Double_t rLow = TMath::Floor(roc->GetPadRowRadii(0,0))-1; // first padRow plus some margin
+
+ // fulcrums in R
+ fgkRList[0] = rLow;
+ for (Int_t i = 1; i245)
+ fgkRList[i] = fgkRList[i-1] + 0.5; // 0.5 cm spacing
+ else if (fgkRList[i]<100 || fgkRList[i]>235)
+ fgkRList[i] = fgkRList[i-1] + 1.5; // 1.5 cm spacing
+ else if (fgkRList[i]<120 || fgkRList[i]>225)
+ fgkRList[i] = fgkRList[i-1] + 2.5; // 2.5 cm spacing
+ }
+
+ // fulcrums in Z
+ fgkZList[0] = -249.5;
+ fgkZList[kNZ-1] = 249.5;
+ for (Int_t j = 1; j248)
+ fgkZList[j] = fgkZList[j-1] + 0.5; // 0.5 cm spacing
+ else if (TMath::Abs(fgkZList[j])<10 || TMath::Abs(fgkZList[j])>235)
+ fgkZList[j] = fgkZList[j-1] + 1.5; // 1.5 cm spacing
+ else if (TMath::Abs(fgkZList[j])<25 || TMath::Abs(fgkZList[j])>225)
+ fgkZList[j] = fgkZList[j-1] + 2.5; // 2.5 cm spacing
+ else
+ fgkZList[j] = fgkZList[j-1] + 4; // 4 cm spacing
+
+ fgkZList[kNZ-j-1] = -fgkZList[j];
+ }
+
+ // fulcrums in phi
+ for (Int_t k = 0; k coef1(rows) ;
+ std::vector coef2(rows) ;
+
+ for ( Int_t i = iOne ; i < rows-1 ; i+=iOne ) {
+ Float_t radius = fgkIFCRadius + i*gridSizeR ;
+ coef1[i] = 1.0 + tempGridSizeR/(2*radius);
+ coef2[i] = 1.0 - tempGridSizeR/(2*radius);
+ }
+
+ TMatrixD sumChargeDensity(rows,columns) ;
+
+ for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
+ Float_t radius = fgkIFCRadius + iOne*gridSizeR ;
+ for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
+ if ( iOne == 1 && jOne == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
+ else {
+ // Add up all enclosed charge density contributions within 1/2 unit in all directions
+ Float_t weight = 0.0 ;
+ Float_t sum = 0.0 ;
+ sumChargeDensity(i,j) = 0.0 ;
+ for ( Int_t ii = i-iOne/2 ; ii <= i+iOne/2 ; ii++ ) {
+ for ( Int_t jj = j-jOne/2 ; jj <= j+jOne/2 ; jj++ ) {
+ if ( ii == i-iOne/2 || ii == i+iOne/2 || jj == j-jOne/2 || jj == j+jOne/2 ) weight = 0.5 ;
+ else
+ weight = 1.0 ;
+ // Note that this is cylindrical geometry
+ sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
+ sum += weight*radius ;
+ }
+ }
+ sumChargeDensity(i,j) /= sum ;
+ }
+ sumChargeDensity(i,j) *= tempGridSizeR*tempGridSizeR; // just saving a step later on
+ }
+ }
+
+ for ( Int_t k = 1 ; k <= iterations; k++ ) {
+ // Solve Poisson's Equation
+ // Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1
+ // as interations increase.
+ Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
+ Float_t overRelaxM1 = overRelax - 1.0 ;
+ Float_t overRelaxtempFourth, overRelaxcoef5 ;
+ overRelaxtempFourth = overRelax * tempFourth ;
+ overRelaxcoef5 = overRelaxM1 / overRelaxtempFourth ;
+
+ for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
+ for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
+
+ arrayV(i,j) = ( coef2[i] * arrayV(i-iOne,j)
+ + tempRatio * ( arrayV(i,j-jOne) + arrayV(i,j+jOne) )
+ - overRelaxcoef5 * arrayV(i,j)
+ + coef1[i] * arrayV(i+iOne,j)
+ + sumChargeDensity(i,j)
+ ) * overRelaxtempFourth;
+ }
+ }
+
+ if ( k == iterations ) {
+ // After full solution is achieved, copy low resolution solution into higher res array
+ for ( Int_t i = iOne ; i < rows-1 ; i += iOne ) {
+ for ( Int_t j = jOne ; j < columns-1 ; j += jOne ) {
+
+ if ( iOne > 1 ) {
+ arrayV(i+iOne/2,j) = ( arrayV(i+iOne,j) + arrayV(i,j) ) / 2 ;
+ if ( i == iOne ) arrayV(i-iOne/2,j) = ( arrayV(0,j) + arrayV(iOne,j) ) / 2 ;
+ }
+ if ( jOne > 1 ) {
+ arrayV(i,j+jOne/2) = ( arrayV(i,j+jOne) + arrayV(i,j) ) / 2 ;
+ if ( j == jOne ) arrayV(i,j-jOne/2) = ( arrayV(i,0) + arrayV(i,jOne) ) / 2 ;
+ }
+ if ( iOne > 1 && jOne > 1 ) {
+ arrayV(i+iOne/2,j+jOne/2) = ( arrayV(i+iOne,j+jOne) + arrayV(i,j) ) / 2 ;
+ if ( i == iOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(0,j-jOne) + arrayV(iOne,j) ) / 2 ;
+ if ( j == jOne ) arrayV(i-iOne/2,j-jOne/2) = ( arrayV(i-iOne,0) + arrayV(i,jOne) ) / 2 ;
+ // Note that this leaves a point at the upper left and lower right corners uninitialized.
+ // -> Not a big deal.
+ }
+
+ }
+ }
+ }
+
+ }
+
+ iOne = iOne / 2 ; if ( iOne < 1 ) iOne = 1 ;
+ jOne = jOne / 2 ; if ( jOne < 1 ) jOne = 1 ;
+
+ sumChargeDensity.Clear();
+ }
+
+ // Differentiate V(r) and solve for E(r) using special equations for the first and last rows
+ for ( Int_t j = 0 ; j < columns ; j++ ) {
+ for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayEr(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
+ arrayEr(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
+ arrayEr(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
+ }
+
+ // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
+ for ( Int_t i = 0 ; i < rows ; i++) {
+ for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayEz(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
+ arrayEz(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
+ arrayEz(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
+ }
+
+ for ( Int_t i = 0 ; i < rows ; i++) {
+ // Note: go back and compare to old version of this code. See notes below.
+ // JT Test ... attempt to divide by real Ez not Ez to first order
+ for ( Int_t j = 0 ; j < columns ; j++ ) {
+ arrayEz(i,j) += ezField;
+ // This adds back the overall Z gradient of the field (main E field component)
+ }
+ // Warning: (-=) assumes you are using an error potetial without the overall Field included
+ }
+
+ // Integrate Er/Ez from Z to zero
+ for ( Int_t j = 0 ; j < columns ; j++ ) {
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+
+ Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
+ arrayErOverEz(i,j) = 0.0 ;
+ arrayDeltaEz(i,j) = 0.0 ;
+
+ for ( Int_t k = j ; k < columns ; k++ ) {
+ arrayErOverEz(i,j) += index*(gridSizeZ/3.0)*arrayEr(i,k)/arrayEz(i,k) ;
+ arrayDeltaEz(i,j) += index*(gridSizeZ/3.0)*(arrayEz(i,k)-ezField) ;
+ if ( index != 4 ) index = 4; else index = 2 ;
+ }
+ if ( index == 4 ) {
+ arrayErOverEz(i,j) -= (gridSizeZ/3.0)*arrayEr(i,columns-1)/arrayEz(i,columns-1) ;
+ arrayDeltaEz(i,j) -= (gridSizeZ/3.0)*(arrayEz(i,columns-1)-ezField) ;
+ }
+ if ( index == 2 ) {
+ arrayErOverEz(i,j) += (gridSizeZ/3.0) * ( 0.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
+ -2.5*arrayEr(i,columns-1)/arrayEz(i,columns-1));
+ arrayDeltaEz(i,j) += (gridSizeZ/3.0) * ( 0.5*(arrayEz(i,columns-2)-ezField)
+ -2.5*(arrayEz(i,columns-1)-ezField));
+ }
+ if ( j == columns-2 ) {
+ arrayErOverEz(i,j) = (gridSizeZ/3.0) * ( 1.5*arrayEr(i,columns-2)/arrayEz(i,columns-2)
+ +1.5*arrayEr(i,columns-1)/arrayEz(i,columns-1) ) ;
+ arrayDeltaEz(i,j) = (gridSizeZ/3.0) * ( 1.5*(arrayEz(i,columns-2)-ezField)
+ +1.5*(arrayEz(i,columns-1)-ezField) ) ;
+ }
+ if ( j == columns-1 ) {
+ arrayErOverEz(i,j) = 0.0 ;
+ arrayDeltaEz(i,j) = 0.0 ;
+ }
+ }
+ }
+
+ // calculate z distortion from the integrated Delta Ez residuals
+ // and include the aquivalence (Volt to cm) of the ROC shift !!
+
+ for ( Int_t j = 0 ; j < columns ; j++ ) {
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+
+ // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
+ arrayDeltaEz(i,j) = arrayDeltaEz(i,j)*fgkdvdE;
+
+ // ROC Potential in cm aquivalent
+ Double_t dzROCShift = arrayV(i, columns -1)/ezField;
+ if ( rocDisplacement ) arrayDeltaEz(i,j) = arrayDeltaEz(i,j) + dzROCShift; // add the ROC misaligment
+
+ }
+ }
+
+ arrayEr.Clear();
+ arrayEz.Clear();
+
+}
+
+void AliTPCCorrection::PoissonRelaxation3D( TMatrixD**arrayofArrayV, TMatrixD**arrayofChargeDensities,
+ TMatrixD**arrayofEroverEz, TMatrixD**arrayofEPhioverEz, TMatrixD**arrayofDeltaEz,
+ Int_t rows, Int_t columns, Int_t phislices,
+ Float_t deltaphi, Int_t iterations, Int_t symmetry,
+ Bool_t rocDisplacement ) {
+ //
+ // 3D - Solve Poisson's Equation in 3D by Relaxation Technique
+ //
+ // NOTE: In order for this algorith to work, the number of rows and columns must be a power of 2 plus one.
+ // The number of rows and COLUMNS can be different.
+ //
+ // ROWS == 2**M + 1
+ // COLUMNS == 2**N + 1
+ // PHISLICES == Arbitrary but greater than 3
+ //
+ // DeltaPhi in Radians
+ //
+ // SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions
+ // = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).
+ //
+ // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
+
+ const Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
+
+ const Float_t gridSizeR = (fgkOFCRadius-fgkIFCRadius) / (rows-1) ;
+ const Float_t gridSizePhi = deltaphi ;
+ const Float_t gridSizeZ = fgkTPCZ0 / (columns-1) ;
+ const Float_t ratioPhi = gridSizeR*gridSizeR / (gridSizePhi*gridSizePhi) ;
+ const Float_t ratioZ = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
+
+ TMatrixD arrayE(rows,columns) ;
+
+ // Check that the number of rows and columns is suitable for a binary expansion
+ if ( !IsPowerOfTwo((rows-1)) ) {
+ AliError("Poisson3DRelaxation - Error in the number of rows. Must be 2**M - 1");
+ return; }
+ if ( !IsPowerOfTwo((columns-1)) ) {
+ AliError("Poisson3DRelaxation - Error in the number of columns. Must be 2**N - 1");
+ return; }
+ if ( phislices <= 3 ) {
+ AliError("Poisson3DRelaxation - Error in the number of phislices. Must be larger than 3");
+ return; }
+ if ( phislices > 1000 ) {
+ AliError("Poisson3D phislices > 1000 is not allowed (nor wise) ");
+ return; }
+
+ // Solve Poisson's equation in cylindrical coordinates by relaxation technique
+ // Allow for different size grid spacing in R and Z directions
+ // Use a binary expansion of the matrix to speed up the solution of the problem
+
+ Int_t loops, mplus, mminus, signplus, signminus ;
+ Int_t ione = (rows-1)/4 ;
+ Int_t jone = (columns-1)/4 ;
+ loops = TMath::Max(ione, jone) ; // Calculate the number of loops for the binary expansion
+ loops = 1 + (int) ( 0.5 + TMath::Log2((double)loops) ) ; // Solve for N in 2**N
+
+ TMatrixD* arrayofSumChargeDensities[1000] ; // Create temporary arrays to store low resolution charge arrays
+
+ for ( Int_t i = 0 ; i < phislices ; i++ ) { arrayofSumChargeDensities[i] = new TMatrixD(rows,columns) ; }
+ AliSysInfo::AddStamp("3DInit", 10,0,0);
+
+ for ( Int_t count = 0 ; count < loops ; count++ ) { // START the master loop and do the binary expansion
+ AliSysInfo::AddStamp("3Diter", 20,count,0);
+
+ Float_t tempgridSizeR = gridSizeR * ione ;
+ Float_t tempratioPhi = ratioPhi * ione * ione ; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
+ Float_t tempratioZ = ratioZ * ione * ione / ( jone * jone ) ;
+
+ std::vector coef1(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+ std::vector coef2(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+ std::vector coef3(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+ std::vector coef4(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+
+ for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
+ Float_t radius = fgkIFCRadius + i*gridSizeR ;
+ coef1[i] = 1.0 + tempgridSizeR/(2*radius);
+ coef2[i] = 1.0 - tempgridSizeR/(2*radius);
+ coef3[i] = tempratioPhi/(radius*radius);
+ coef4[i] = 0.5 / (1.0 + tempratioZ + coef3[i]);
+ }
+
+ for ( Int_t m = 0 ; m < phislices ; m++ ) {
+ TMatrixD &chargeDensity = *arrayofChargeDensities[m] ;
+ TMatrixD &sumChargeDensity = *arrayofSumChargeDensities[m] ;
+ for ( Int_t i = ione ; i < rows-1 ; i += ione ) {
+ Float_t radius = fgkIFCRadius + i*gridSizeR ;
+ for ( Int_t j = jone ; j < columns-1 ; j += jone ) {
+ if ( ione == 1 && jone == 1 ) sumChargeDensity(i,j) = chargeDensity(i,j) ;
+ else { // Add up all enclosed charge density contributions within 1/2 unit in all directions
+ Float_t weight = 0.0 ;
+ Float_t sum = 0.0 ;
+ sumChargeDensity(i,j) = 0.0 ;
+ for ( Int_t ii = i-ione/2 ; ii <= i+ione/2 ; ii++ ) {
+ for ( Int_t jj = j-jone/2 ; jj <= j+jone/2 ; jj++ ) {
+ if ( ii == i-ione/2 || ii == i+ione/2 || jj == j-jone/2 || jj == j+jone/2 ) weight = 0.5 ;
+ else
+ weight = 1.0 ;
+ sumChargeDensity(i,j) += chargeDensity(ii,jj)*weight*radius ;
+ sum += weight*radius ;
+ }
+ }
+ sumChargeDensity(i,j) /= sum ;
+ }
+ sumChargeDensity(i,j) *= tempgridSizeR*tempgridSizeR; // just saving a step later on
+ }
+ }
+ }
+
+ for ( Int_t k = 1 ; k <= iterations; k++ ) {
+
+ // over-relaxation index, >= 1 but < 2
+ Float_t overRelax = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/iterations ) ) ;
+ Float_t overRelaxM1 = overRelax - 1.0 ;
+
+ std::vector overRelaxcoef4(rows) ; // Do this the standard C++ way to avoid gcc extensions
+ std::vector overRelaxcoef5(rows) ; // Do this the standard C++ way to avoid gcc extensions
+
+ for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
+ overRelaxcoef4[i] = overRelax * coef4[i] ;
+ overRelaxcoef5[i] = overRelaxM1 / overRelaxcoef4[i] ;
+ }
+
+ for ( Int_t m = 0 ; m < phislices ; m++ ) {
+
+ mplus = m + 1; signplus = 1 ;
+ mminus = m - 1 ; signminus = 1 ;
+ if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
+ if ( mplus > phislices-1 ) mplus = phislices - 2 ;
+ if ( mminus < 0 ) mminus = 1 ;
+ }
+ else if (symmetry==-1) { // Anti-symmetry in phi
+ if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
+ if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
+ }
+ else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
+ if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
+ if ( mminus < 0 ) mminus = m - 1 + phislices ;
+ }
+ TMatrixD& arrayV = *arrayofArrayV[m] ;
+ TMatrixD& arrayVP = *arrayofArrayV[mplus] ;
+ TMatrixD& arrayVM = *arrayofArrayV[mminus] ;
+ TMatrixD& sumChargeDensity = *arrayofSumChargeDensities[m] ;
+ Double_t *arrayVfast = arrayV.GetMatrixArray();
+ Double_t *arrayVPfast = arrayVP.GetMatrixArray();
+ Double_t *arrayVMfast = arrayVM.GetMatrixArray();
+ Double_t *sumChargeDensityFast=sumChargeDensity.GetMatrixArray();
+
+ if (0){
+ // slow implementation
+ for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
+ for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
+
+ arrayV(i,j) = ( coef2[i] * arrayV(i-ione,j)
+ + tempratioZ * ( arrayV(i,j-jone) + arrayV(i,j+jone) )
+ - overRelaxcoef5[i] * arrayV(i,j)
+ + coef1[i] * arrayV(i+ione,j)
+ + coef3[i] * ( signplus*arrayVP(i,j) + signminus*arrayVM(i,j) )
+ + sumChargeDensity(i,j)
+ ) * overRelaxcoef4[i] ;
+ // Note: over-relax the solution at each step. This speeds up the convergance.
+ }
+ }
+ }else{
+ for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
+ Double_t *arrayVfastI = &(arrayVfast[i*columns]);
+ Double_t *arrayVPfastI = &(arrayVPfast[i*columns]);
+ Double_t *arrayVMfastI = &(arrayVMfast[i*columns]);
+ Double_t *sumChargeDensityFastI=&(sumChargeDensityFast[i*columns]);
+ for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
+ Double_t /*resSlow*/resFast;
+// resSlow = ( coef2[i] * arrayV(i-ione,j)
+// + tempratioZ * ( arrayV(i,j-jone) + arrayV(i,j+jone) )
+// - overRelaxcoef5[i] * arrayV(i,j)
+// + coef1[i] * arrayV(i+ione,j)
+// + coef3[i] * ( signplus*arrayVP(i,j) + signminus*arrayVM(i,j) )
+// + sumChargeDensity(i,j)
+// ) * overRelaxcoef4[i] ;
+ resFast = ( coef2[i] * arrayVfastI[j-columns*ione]
+ + tempratioZ * ( arrayVfastI[j-jone] + arrayVfastI[j+jone] )
+ - overRelaxcoef5[i] * arrayVfastI[j]
+ + coef1[i] * arrayVfastI[j+columns*ione]
+ + coef3[i] * ( signplus* arrayVPfastI[j] + signminus*arrayVMfastI[j])
+ + sumChargeDensityFastI[j]
+ ) * overRelaxcoef4[i] ;
+// if (resSlow!=resFast){
+// printf("problem\t%d\t%d\t%f\t%f\t%f\n",i,j,resFast,resSlow,resFast-resSlow);
+// }
+ arrayVfastI[j]=resFast;
+ // Note: over-relax the solution at each step. This speeds up the convergance.
+ }
+ }
+ }
+
+ if ( k == iterations ) { // After full solution is achieved, copy low resolution solution into higher res array
+ for ( Int_t i = ione ; i < rows-1 ; i+=ione ) {
+ for ( Int_t j = jone ; j < columns-1 ; j+=jone ) {
+
+ if ( ione > 1 ) {
+ arrayV(i+ione/2,j) = ( arrayV(i+ione,j) + arrayV(i,j) ) / 2 ;
+ if ( i == ione ) arrayV(i-ione/2,j) = ( arrayV(0,j) + arrayV(ione,j) ) / 2 ;
+ }
+ if ( jone > 1 ) {
+ arrayV(i,j+jone/2) = ( arrayV(i,j+jone) + arrayV(i,j) ) / 2 ;
+ if ( j == jone ) arrayV(i,j-jone/2) = ( arrayV(i,0) + arrayV(i,jone) ) / 2 ;
+ }
+ if ( ione > 1 && jone > 1 ) {
+ arrayV(i+ione/2,j+jone/2) = ( arrayV(i+ione,j+jone) + arrayV(i,j) ) / 2 ;
+ if ( i == ione ) arrayV(i-ione/2,j-jone/2) = ( arrayV(0,j-jone) + arrayV(ione,j) ) / 2 ;
+ if ( j == jone ) arrayV(i-ione/2,j-jone/2) = ( arrayV(i-ione,0) + arrayV(i,jone) ) / 2 ;
+ // Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.
+ }
+ }
+ }
+ }
+
+ }
+ }
+
+ ione = ione / 2 ; if ( ione < 1 ) ione = 1 ;
+ jone = jone / 2 ; if ( jone < 1 ) jone = 1 ;
+
+ }
+
+ //Differentiate V(r) and solve for E(r) using special equations for the first and last row
+ //Integrate E(r)/E(z) from point of origin to pad plane
+ AliSysInfo::AddStamp("CalcField", 100,0,0);
+
+ for ( Int_t m = 0 ; m < phislices ; m++ ) {
+ TMatrixD& arrayV = *arrayofArrayV[m] ;
+ TMatrixD& eroverEz = *arrayofEroverEz[m] ;
+
+ for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
+
+ // Differentiate in R
+ for ( Int_t i = 1 ; i < rows-1 ; i++ ) arrayE(i,j) = -1 * ( arrayV(i+1,j) - arrayV(i-1,j) ) / (2*gridSizeR) ;
+ arrayE(0,j) = -1 * ( -0.5*arrayV(2,j) + 2.0*arrayV(1,j) - 1.5*arrayV(0,j) ) / gridSizeR ;
+ arrayE(rows-1,j) = -1 * ( 1.5*arrayV(rows-1,j) - 2.0*arrayV(rows-2,j) + 0.5*arrayV(rows-3,j) ) / gridSizeR ;
+ // Integrate over Z
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+ Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
+ eroverEz(i,j) = 0.0 ;
+ for ( Int_t k = j ; k < columns ; k++ ) {
+
+ eroverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
+ if ( index != 4 ) index = 4; else index = 2 ;
+ }
+ if ( index == 4 ) eroverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
+ if ( index == 2 ) eroverEz(i,j) +=
+ (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
+ if ( j == columns-2 ) eroverEz(i,j) =
+ (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
+ if ( j == columns-1 ) eroverEz(i,j) = 0.0 ;
+ }
+ }
+ // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
+ // eroverEz.Draw("surf") ; } // JT test
+ }
+ AliSysInfo::AddStamp("IntegrateEr", 120,0,0);
+
+ //Differentiate V(r) and solve for E(phi)
+ //Integrate E(phi)/E(z) from point of origin to pad plane
+
+ for ( Int_t m = 0 ; m < phislices ; m++ ) {
+
+ mplus = m + 1; signplus = 1 ;
+ mminus = m - 1 ; signminus = 1 ;
+ if (symmetry==1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
+ if ( mplus > phislices-1 ) mplus = phislices - 2 ;
+ if ( mminus < 0 ) mminus = 1 ;
+ }
+ else if (symmetry==-1) { // Anti-symmetry in phi
+ if ( mplus > phislices-1 ) { mplus = phislices - 2 ; signplus = -1 ; }
+ if ( mminus < 0 ) { mminus = 1 ; signminus = -1 ; }
+ }
+ else { // No Symmetries in phi, no boundaries, the calculations is continuous across all phi
+ if ( mplus > phislices-1 ) mplus = m + 1 - phislices ;
+ if ( mminus < 0 ) mminus = m - 1 + phislices ;
+ }
+ TMatrixD &arrayVP = *arrayofArrayV[mplus] ;
+ TMatrixD &arrayVM = *arrayofArrayV[mminus] ;
+ TMatrixD &ePhioverEz = *arrayofEPhioverEz[m] ;
+ for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
+ // Differentiate in Phi
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+ Float_t radius = fgkIFCRadius + i*gridSizeR ;
+ arrayE(i,j) = -1 * (signplus * arrayVP(i,j) - signminus * arrayVM(i,j) ) / (2*radius*gridSizePhi) ;
+ }
+ // Integrate over Z
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+ Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
+ ePhioverEz(i,j) = 0.0 ;
+ for ( Int_t k = j ; k < columns ; k++ ) {
+
+ ePhioverEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k)/(-1*ezField) ;
+ if ( index != 4 ) index = 4; else index = 2 ;
+ }
+ if ( index == 4 ) ePhioverEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1)/ (-1*ezField) ;
+ if ( index == 2 ) ePhioverEz(i,j) +=
+ (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1))/(-1*ezField) ;
+ if ( j == columns-2 ) ePhioverEz(i,j) =
+ (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1))/(-1*ezField) ;
+ if ( j == columns-1 ) ePhioverEz(i,j) = 0.0 ;
+ }
+ }
+ // if ( m == 5 ) { TCanvas* c2 = new TCanvas("arrayE","arrayE",50,50,840,600) ; c2 -> cd() ;
+ // arrayE.Draw("surf") ; } // JT test
+ }
+ AliSysInfo::AddStamp("IntegrateEphi", 130,0,0);
+
+
+ // Differentiate V(r) and solve for E(z) using special equations for the first and last row
+ // Integrate (E(z)-Ezstd) from point of origin to pad plane
+
+ for ( Int_t m = 0 ; m < phislices ; m++ ) {
+ TMatrixD& arrayV = *arrayofArrayV[m] ;
+ TMatrixD& deltaEz = *arrayofDeltaEz[m] ;
+
+ // Differentiate V(z) and solve for E(z) using special equations for the first and last columns
+ for ( Int_t i = 0 ; i < rows ; i++) {
+ for ( Int_t j = 1 ; j < columns-1 ; j++ ) arrayE(i,j) = -1 * ( arrayV(i,j+1) - arrayV(i,j-1) ) / (2*gridSizeZ) ;
+ arrayE(i,0) = -1 * ( -0.5*arrayV(i,2) + 2.0*arrayV(i,1) - 1.5*arrayV(i,0) ) / gridSizeZ ;
+ arrayE(i,columns-1) = -1 * ( 1.5*arrayV(i,columns-1) - 2.0*arrayV(i,columns-2) + 0.5*arrayV(i,columns-3) ) / gridSizeZ ;
+ }
+
+ for ( Int_t j = columns-1 ; j >= 0 ; j-- ) { // Count backwards to facilitate integration over Z
+ // Integrate over Z
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+ Int_t index = 1 ; // Simpsons rule if N=odd. If N!=odd then add extra point by trapezoidal rule.
+ deltaEz(i,j) = 0.0 ;
+ for ( Int_t k = j ; k < columns ; k++ ) {
+ deltaEz(i,j) += index*(gridSizeZ/3.0)*arrayE(i,k) ;
+ if ( index != 4 ) index = 4; else index = 2 ;
+ }
+ if ( index == 4 ) deltaEz(i,j) -= (gridSizeZ/3.0)*arrayE(i,columns-1) ;
+ if ( index == 2 ) deltaEz(i,j) +=
+ (gridSizeZ/3.0)*(0.5*arrayE(i,columns-2)-2.5*arrayE(i,columns-1)) ;
+ if ( j == columns-2 ) deltaEz(i,j) =
+ (gridSizeZ/3.0)*(1.5*arrayE(i,columns-2)+1.5*arrayE(i,columns-1)) ;
+ if ( j == columns-1 ) deltaEz(i,j) = 0.0 ;
+ }
+ }
+
+ // if ( m == 0 ) { TCanvas* c1 = new TCanvas("erOverEz","erOverEz",50,50,840,600) ; c1 -> cd() ;
+ // eroverEz.Draw("surf") ; } // JT test
+
+ // calculate z distortion from the integrated Delta Ez residuals
+ // and include the aquivalence (Volt to cm) of the ROC shift !!
+
+ for ( Int_t j = 0 ; j < columns ; j++ ) {
+ for ( Int_t i = 0 ; i < rows ; i++ ) {
+
+ // Scale the Ez distortions with the drift velocity pertubation -> delivers cm
+ deltaEz(i,j) = deltaEz(i,j)*fgkdvdE;
+
+ // ROC Potential in cm aquivalent
+ Double_t dzROCShift = arrayV(i, columns -1)/ezField;
+ if ( rocDisplacement ) deltaEz(i,j) = deltaEz(i,j) + dzROCShift; // add the ROC misaligment
+
+ }
+ }
+
+ } // end loop over phi
+ AliSysInfo::AddStamp("IntegrateEz", 140,0,0);
+
+
+ for ( Int_t k = 0 ; k < phislices ; k++ )
+ {
+ arrayofSumChargeDensities[k]->Delete() ;
+ }
+
+
+
+ arrayE.Clear();
+}
+
+
+Int_t AliTPCCorrection::IsPowerOfTwo(Int_t i) const {
+ //
+ // Helperfunction: Check if integer is a power of 2
+ //
+ Int_t j = 0;
+ while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
+ if ( j == 1 ) return(1) ; // True
+ return(0) ; // False
+}
+
+
+AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
+ //
+ // Fit the track parameters - without and with distortion
+ // 1. Space points in the TPC are simulated along the trajectory
+ // 2. Space points distorted
+ // 3. Fits the non distorted and distroted track to the reference plane at refX
+ // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
+ //
+ // trackIn - input track parameters
+ // refX - reference X to fit the track
+ // dir - direction - out=1 or in=-1
+ // pcstream - debug streamer to check the results
+ //
+ // see AliExternalTrackParam.h documentation:
+ // track1.fP[0] - local y (rphi)
+ // track1.fP[1] - z
+ // track1.fP[2] - sinus of local inclination angle
+ // track1.fP[3] - tangent of deep angle
+ // track1.fP[4] - 1/pt
+
+ AliTPCROC * roc = AliTPCROC::Instance();
+ const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
+ const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
+ const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
+ const Double_t kMaxSnp = 0.85;
+ const Double_t kSigmaY=0.1;
+ const Double_t kSigmaZ=0.1;
+ const Double_t kMaxR=500;
+ const Double_t kMaxZ=500;
+
+ const Double_t kMaxZ0=220;
+ const Double_t kZcut=3;
+ const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+ Int_t npoints1=0;
+ Int_t npoints2=0;
+
+ AliExternalTrackParam track(trackIn); //
+ // generate points
+ AliTrackPointArray pointArray0(npoints0);
+ AliTrackPointArray pointArray1(npoints0);
+ Double_t xyz[3];
+ if (!AliTrackerBase::PropagateTrackTo(&track,kRTPC0,kMass,5,kTRUE,kMaxSnp)) return 0;
+ //
+ // simulate the track
+ Int_t npoints=0;
+ Float_t covPoint[6]={0,0,0, static_cast(kSigmaY*kSigmaY),0,static_cast(kSigmaZ*kSigmaZ)}; //covariance at the local frame
+ for (Double_t radius=kRTPC0; radiusGaus(0,0.000005);
+ xyz[1]+=gRandom->Gaus(0,0.000005);
+ xyz[2]+=gRandom->Gaus(0,0.000005);
+ if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
+ if (TMath::Abs(track.GetX())kRTPC1) continue;
+ AliTrackPoint pIn0; // space point
+ AliTrackPoint pIn1;
+ Int_t sector= (xyz[2]>0)? 0:18;
+ pointArray0.GetPoint(pIn0,npoints);
+ pointArray1.GetPoint(pIn1,npoints);
+ Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
+ Float_t distPoint[3]={static_cast(xyz[0]),static_cast(xyz[1]),static_cast(xyz[2])};
+ DistortPoint(distPoint, sector);
+ pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
+ pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
+ //
+ track.Rotate(alpha);
+ AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
+ AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
+ prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
+ prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
+ pIn0=prot0.Rotate(-alpha); // rotate back to global frame
+ pIn1=prot1.Rotate(-alpha); // rotate back to global frame
+ pointArray0.AddPoint(npoints, &pIn0);
+ pointArray1.AddPoint(npoints, &pIn1);
+ npoints++;
+ if (npoints>=npoints0) break;
+ }
+ if (npointsResetCovariance(10);
+ track1->ResetCovariance(10);
+ if (TMath::Abs(AliTrackerBase::GetBz())<0.01){
+ ((Double_t*)track0->GetParameter())[4]= trackIn.GetParameter()[4];
+ ((Double_t*)track1->GetParameter())[4]= trackIn.GetParameter()[4];
+ }
+ for (Int_t jpoint=0; jpoint0) ? jpoint: npoints-1-jpoint;
+ //
+ AliTrackPoint pIn0;
+ AliTrackPoint pIn1;
+ pointArray0.GetPoint(pIn0,ipoint);
+ pointArray1.GetPoint(pIn1,ipoint);
+ AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
+ AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
+ if (TMath::Abs(prot0.GetX())kRTPC1) continue;
+ //
+ if (!AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
+ if (!AliTrackerBase::PropagateTrackTo(track1,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
+ if (TMath::Abs(track0->GetZ())>kMaxZ) break;
+ if (TMath::Abs(track0->GetX())>kMaxR) break;
+ if (TMath::Abs(track1->GetZ())>kMaxZ) break;
+ if (TMath::Abs(track1->GetX())>kMaxR) break;
+ if (dir>0 && track1->GetX()>refX) continue;
+ if (dir<0 && track1->GetX()GetZ())Update(pointPos,pointCov)) break;
+ //
+ Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
+ Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
+ Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
+
+ pointPos[0]=prot1.GetY()-deltaYX;//local y is sign correct? should be minus
+ pointPos[1]=prot1.GetZ()-deltaZX;//local z is sign correct? should be minus
+ pointCov[0]=prot1.GetCov()[3];//simay^2
+ pointCov[1]=prot1.GetCov()[4];//sigmayz
+ pointCov[2]=prot1.GetCov()[5];//sigmaz^2
+ if (!track1->Update(pointPos,pointCov)) break;
+ npoints1++;
+ npoints2++;
+ }
+ if (npoints2Rotate(track0->GetAlpha());
+ AliTrackerBase::PropagateTrackTo(track1,track0->GetX(),kMass,5.,kFALSE,kMaxSnp);
+
+ if (pcstream) (*pcstream)<