+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; i<kNR; i++) {
+ fgkRList[i] = fgkRList[i-1] + 3.5; // 3.5 cm spacing
+ if (fgkRList[i]<90 ||fgkRList[i]>245)
+ 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; j<kNZ/2; j++) {
+ fgkZList[j] = fgkZList[j-1];
+ if (TMath::Abs(fgkZList[j])< 0.15)
+ fgkZList[j] = fgkZList[j-1] + 0.09; // 0.09 cm spacing
+ else if(TMath::Abs(fgkZList[j])< 0.6)
+ fgkZList[j] = fgkZList[j-1] + 0.4; // 0.4 cm spacing
+ else if (TMath::Abs(fgkZList[j])< 2.5 || TMath::Abs(fgkZList[j])>248)
+ 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<kNPhi; k++)
+ fgkPhiList[k] = TMath::TwoPi()*k/(kNPhi-1);
+
+
+}
+
+
+void AliTPCCorrection::PoissonRelaxation2D(TMatrixD &arrayV, TMatrixD &chargeDensity,
+ TMatrixD &arrayErOverEz, TMatrixD &arrayDeltaEz,
+ const Int_t rows, const Int_t columns, const Int_t iterations,
+ const Bool_t rocDisplacement ) {
+ //
+ // Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)
+ //
+ // Solve Poissons equation in a cylindrical coordinate system. The arrayV matrix must be filled with the
+ // boundary conditions on the first and last rows, and the first and last columns. The remainder of the
+ // array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains
+ // the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if
+ // you wish to solve Laplaces equation however it should not contain random numbers or you will get
+ // random numbers back as a solution.
+ // Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to
+ // speed up the convergence to the best solution, this algorithm does a binary expansion of the solution
+ // space. First it solves the problem on a very sparse grid by skipping rows and columns in the original
+ // matrix. Then it doubles the number of points and solves the problem again. Then it doubles the
+ // number of points and solves the problem again. This happens several times until the maximum number
+ // of points has been included in the array.
+ //
+ // NOTE: In order for this algorithmto work, the number of rows and columns must be a power of 2 plus one.
+ // So rows == 2**M + 1 and columns == 2**N + 1. The number of rows and columns can be different.
+ //
+ // NOTE: rocDisplacement is used to include (or ignore) the ROC misalignment in the dz calculation
+ //
+ // Original code by Jim Thomas (STAR TPC Collaboration)
+ //
+
+ 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 gridSizeZ = fgkTPCZ0 / (columns-1) ;
+ const Float_t ratio = gridSizeR*gridSizeR / (gridSizeZ*gridSizeZ) ;
+
+ TMatrixD arrayEr(rows,columns) ;
+ TMatrixD arrayEz(rows,columns) ;
+
+ //Check that number of rows and columns is suitable for a binary expansion
+
+ if ( !IsPowerOfTwo(rows-1) ) {
+ AliError("PoissonRelaxation - Error in the number of rows. Must be 2**M - 1");
+ return;
+ }
+ if ( !IsPowerOfTwo(columns-1) ) {
+ AliError("PoissonRelaxation - Error in the number of columns. Must be 2**N - 1");
+ 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 size of the matrix to speed up the solution of the problem
+
+ Int_t iOne = (rows-1)/4 ;
+ Int_t jOne = (columns-1)/4 ;
+ // Solve for N in 2**N, add one.
+ Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (double) TMath::Max(iOne,jOne) ) ) ;
+
+ for ( Int_t count = 0 ; count < loops ; count++ ) {
+ // Loop while the matrix expands & the resolution increases.
+
+ Float_t tempGridSizeR = gridSizeR * iOne ;
+ Float_t tempRatio = ratio * iOne * iOne / ( jOne * jOne ) ;
+ Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
+
+ // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+ std::vector<float> coef1(rows) ;
+ std::vector<float> 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,
+ const Int_t rows, const Int_t columns, const Int_t phislices,
+ const Float_t deltaphi, const Int_t iterations, const 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) ; }
+
+ for ( Int_t count = 0 ; count < loops ; count++ ) { // START the master loop and do the binary expansion
+
+ 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<float> coef1(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+ std::vector<float> coef2(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+ std::vector<float> coef3(rows) ; // Do this the standard C++ way to avoid gcc extensions for Float_t coef1[rows]
+ std::vector<float> 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<float> overRelaxcoef4(rows) ; // Do this the standard C++ way to avoid gcc extensions
+ std::vector<float> 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] ;
+
+ 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.
+
+ }
+ }
+
+ 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
+
+ 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
+ }
+
+ //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
+ }
+
+
+ // 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
+
+
+
+ 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
+}
+