Fixes for coding violations
[u/mrichter/AliRoot.git] / TPC / AliTPCFCVoltError3D.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // _________________________________________________________________
17 //
18 // Begin_Html
19 //   <h2>AliTPCFCVoltError3D class   </h2>       
20 //   The class calculates the space point distortions due to residual voltage errors 
21 //   on the Field Cage (FC) boundaries (rods and strips) of the TPC in 3D. It uses 
22 //   the Poisson relaxation technique in three dimension as implemented in the parent class. 
23 //   <p>
24 //   Although the calculation is performed in 3D, the calculation time was significantly 
25 //   reduced by using certain symmetry conditions within the calculation.
26 //   <p>
27 //   The input parameters can be set via the functions (e.g.) SetRodVoltShift(rod,dV) where 
28 //   rod is the number of the rod on which the voltage offset dV is set. The corresponding 
29 //   shift in z direction would be $dz=dV/40$ with an opposite sign for the C side. The 
30 //   rods are numbered in anti-clock-wise direction starting at $\phi=0$. Rods in the IFC 
31 //   are from 0 to 17, rods on the OFC go from 18 to 35. <br>
32 //   This convention is following the offline numbering scheme of the ROCs.
33 //   <p>
34 //   Different misalignment scenarios can be modeled: 
35 //   <ul>
36 //   <li> A rotated clip scenario is only possible at two positions (Rod 11 at IFC, rod 3(+18) at OFC) 
37 //        and can be set via SetRotatedClipVolt. The key feature is that at the mentioned positions,
38 //        the strip-ends were combined. At these positions, a systematic offset of one strip-end in
39 //        respect to the other is possible. 
40 //   <li> A normal rod offset, where the strips follow the movement of the rod, can be set via 
41 //        SetRodVoltShift. It has a anti-mirrored signature in terms of distortions when compared 
42 //        to the rotated clip. This misalignment is possible at each single rod of the FC.
43 //   <li> A simple rod offset, where the strips do not follow the shift, results in an even more 
44 //        localized distortion close to the rod. The difference to a rod shift, where the strips follow,
45 //        is more dominant on the OFC. This effect can be set via the function SetCopperRodShift.
46 //   </ul>
47 // End_Html
48 //
49 // Begin_Macro(source) 
50 //   {
51 //   gROOT->SetStyle("Plain"); gStyle->SetPalette(1);
52 //   TCanvas *c2 = new TCanvas("c2","c2",500,450); 
53 //   AliTPCFCVoltError3D fc;
54 //   fc.SetOmegaTauT1T2(0,1,1); 
55 //   fc.SetRotatedClipVoltA(0,40);
56 //   fc.SetRodVoltShiftA(3,40); 
57 //   fc.SetCopperRodShiftA(7+18,40);
58 //   fc.SetRodVoltShiftA(15+18,40); 
59 //   fc.CreateHistoDRPhiinXY(10)->Draw("cont4z");
60 //   return c2;
61 //   } 
62 // End_Macro
63 //
64 // Begin_Html
65 //   <p>
66 //   Date: 08/08/2010  <br>
67 //   Authors: Jim Thomas, Stefan Rossegger  
68 // End_Html 
69 // _________________________________________________________________
70
71
72 #include "AliMagF.h"
73 #include "TGeoGlobalMagField.h"
74 #include "AliTPCcalibDB.h"
75 #include "AliTPCParam.h"
76 #include "AliLog.h"
77 #include "TMatrixD.h"
78 #include "TMatrixF.h"
79
80 #include "TMath.h"
81 #include "AliTPCROC.h"
82 #include "AliTPCFCVoltError3D.h"
83
84 ClassImp(AliTPCFCVoltError3D)
85
86 AliTPCFCVoltError3D::AliTPCFCVoltError3D()
87   : AliTPCCorrection("FieldCageVoltErrors","FieldCage (Rods) Voltage Errors"),
88     fC0(0.),fC1(0.),
89     fInitLookUp(kFALSE)
90 {
91   //
92   // default constructor
93   //
94
95   // flags for filled 'basic lookup tables'
96   for (Int_t i=0; i<6; i++){
97     fInitLookUpBasic[i]= kFALSE;  
98   }
99
100   // Boundary settings 
101   for (Int_t i=0; i<36; i++){
102     fRodVoltShiftA[i] = 0;  
103     fRodVoltShiftC[i] = 0;  
104   }
105   for (Int_t i=0; i<2; i++){
106     fRotatedClipVoltA[i] = 0;  
107     fRotatedClipVoltC[i] = 0;  
108   }
109   // 
110   for (Int_t i=0; i<36; i++){
111     fCopperRodShiftA[i] = 0;  
112     fCopperRodShiftC[i] = 0;  
113   }
114
115   // Array which will contain the solution according to the setted boundary conditions
116   // it represents a sum up of the 4 basic look up tables (created individually)
117   // see InitFCVoltError3D() function
118   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
119     fLookUpErOverEz[k]   =  new TMatrixF(kNR,kNZ);  
120     fLookUpEphiOverEz[k] =  new TMatrixF(kNR,kNZ);
121     fLookUpDeltaEz[k]    =  new TMatrixF(kNR,kNZ);   
122   }
123   
124   for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
125     fLookUpBasic1ErOverEz[k]   = 0;
126     fLookUpBasic1EphiOverEz[k] = 0; 
127     fLookUpBasic1DeltaEz[k]    = 0;
128
129     fLookUpBasic2ErOverEz[k]   = 0;
130     fLookUpBasic2EphiOverEz[k] = 0; 
131     fLookUpBasic2DeltaEz[k]    = 0;
132
133     fLookUpBasic3ErOverEz[k]   = 0;
134     fLookUpBasic3EphiOverEz[k] = 0; 
135     fLookUpBasic3DeltaEz[k]    = 0;
136
137     fLookUpBasic4ErOverEz[k]   = 0;
138     fLookUpBasic4EphiOverEz[k] = 0; 
139     fLookUpBasic4DeltaEz[k]    = 0;
140     
141     fLookUpBasic5ErOverEz[k]   = 0;
142     fLookUpBasic5EphiOverEz[k] = 0; 
143     fLookUpBasic5DeltaEz[k]    = 0;
144
145     fLookUpBasic6ErOverEz[k]   = 0;
146     fLookUpBasic6EphiOverEz[k] = 0; 
147     fLookUpBasic6DeltaEz[k]    = 0;
148   }
149
150 }
151
152 AliTPCFCVoltError3D::~AliTPCFCVoltError3D() {
153   //
154   // destructor
155   //
156   
157   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
158     delete fLookUpErOverEz[k];
159     delete fLookUpEphiOverEz[k];
160     delete fLookUpDeltaEz[k];
161   }
162
163   for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
164     delete fLookUpBasic1ErOverEz[k];  // does nothing if pointer is zero!
165     delete fLookUpBasic1EphiOverEz[k]; 
166     delete fLookUpBasic1DeltaEz[k]; 
167
168     delete fLookUpBasic2ErOverEz[k];  // does nothing if pointer is zero!
169     delete fLookUpBasic2EphiOverEz[k]; 
170     delete fLookUpBasic2DeltaEz[k]; 
171     
172     delete fLookUpBasic3ErOverEz[k];  // does nothing if pointer is zero!
173     delete fLookUpBasic3EphiOverEz[k]; 
174     delete fLookUpBasic3DeltaEz[k]; 
175
176     delete fLookUpBasic4ErOverEz[k];  // does nothing if pointer is zero!
177     delete fLookUpBasic4EphiOverEz[k]; 
178     delete fLookUpBasic4DeltaEz[k]; 
179
180     delete fLookUpBasic5ErOverEz[k];  // does nothing if pointer is zero!
181     delete fLookUpBasic5EphiOverEz[k]; 
182     delete fLookUpBasic5DeltaEz[k]; 
183
184     delete fLookUpBasic6ErOverEz[k];  // does nothing if pointer is zero!
185     delete fLookUpBasic6EphiOverEz[k]; 
186     delete fLookUpBasic6DeltaEz[k]; 
187
188   }
189 }
190
191 void AliTPCFCVoltError3D::Init() {
192   //
193   // Initialization funtion
194   //
195   
196   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
197   if (!magF) AliError("Magneticd field - not initialized");
198   Double_t bzField = magF->SolenoidField()/10.; //field in T
199   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
200   if (!param) AliError("Parameters - not initialized");
201   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
202   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
203   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 
204   // Correction Terms for effective omegaTau; obtained by a laser calibration run
205   SetOmegaTauT1T2(wt,fT1,fT2);
206
207   if (!fInitLookUp) InitFCVoltError3D();
208 }
209
210 void AliTPCFCVoltError3D::Update(const TTimeStamp &/*timeStamp*/) {
211   //
212   // Update function 
213   //
214   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
215   if (!magF) AliError("Magneticd field - not initialized");
216   Double_t bzField = magF->SolenoidField()/10.; //field in T
217   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
218   if (!param) AliError("Parameters - not initialized");
219   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]  // From dataBase: to be updated: per second (ideally)
220   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
221   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 
222   // Correction Terms for effective omegaTau; obtained by a laser calibration run
223   SetOmegaTauT1T2(wt,fT1,fT2);
224
225
226 }
227
228
229
230 void AliTPCFCVoltError3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
231   //
232   // Calculates the correction due e.g. residual voltage errors on the TPC boundaries
233   //   
234   const Double_t kEpsilon=Double_t(FLT_MIN);
235
236   if (!fInitLookUp) {
237     AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
238     InitFCVoltError3D();
239   }
240
241   static Bool_t forceInit=kTRUE; //temporary needed for back compatibility with old OCDB entries
242   if (forceInit &&fLookUpErOverEz[0]){
243     if (TMath::Abs(fLookUpErOverEz[0]->Sum())<kEpsilon){//temporary needed for back compatibility with old OCDB entries
244       ForceInitFCVoltError3D();
245     }
246     forceInit=kFALSE;
247   }
248
249
250   Int_t   order     = 1 ;               // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2         
251                                         // note that the poisson solution was linearly mirroed on this grid!
252   Float_t intEr, intEphi, intDeltaEz;
253   Float_t r, phi, z ;
254   Int_t    sign;
255
256   r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
257   phi    =  TMath::ATan2(x[1],x[0]) ;
258   if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
259   z      =  x[2] ;                                         // Create temporary copy of x[2]
260
261   if ( (roc%36) < 18 ) {
262     sign =  1;       // (TPC A side)
263   } else {
264     sign = -1;       // (TPC C side)
265   }
266   
267   if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
268   if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
269   
270
271   if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
272     AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
273
274   // Get the Er and Ephi field integrals plus the integral over DeltaEz 
275   intEr      = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi, 
276                                   fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz  );
277   intEphi    = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi, 
278                                   fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz  );
279   intDeltaEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi, 
280                                   fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz  );
281
282   //  printf("%lf %lf %lf\n",intEr,intEphi,intDeltaEz);
283
284   // Calculate distorted position
285   if ( r > 0.0 ) {
286     phi =  phi + ( fC0*intEphi - fC1*intEr ) / r;      
287     r   =  r   + ( fC0*intEr   + fC1*intEphi );  
288   }
289   
290   // Calculate correction in cartesian coordinates
291   dx[0] = r * TMath::Cos(phi) - x[0];
292   dx[1] = r * TMath::Sin(phi) - x[1]; 
293   dx[2] = intDeltaEz;  // z distortion - (internally scaled with driftvelocity dependency 
294                        // on the Ez field plus the actual ROC misalignment (if set TRUE)
295
296 }
297
298 void AliTPCFCVoltError3D::InitFCVoltError3D() {
299   //
300   // Initialization of the Lookup table which contains the solutions of the 
301   // Dirichlet boundary problem
302   // Calculation of the single 3D-Poisson solver is done just if needed
303   // (see basic lookup tables in header file)
304   //
305
306   const Int_t   order       =    1  ;  // Linear interpolation = 1, Quadratic = 2  
307   const Float_t gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
308   const Float_t gridSizeZ   =  fgkTPCZ0 / (kColumns-1) ;
309   const Float_t gridSizePhi =  TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);
310
311   // temporary arrays to create the boundary conditions
312   TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ; 
313   for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
314     arrayofArrayV[k]     =   new TMatrixD(kRows,kColumns) ;
315     arrayofCharge[k]     =   new TMatrixD(kRows,kColumns) ;
316   }
317   Double_t  innerList[kPhiSlices], outerList[kPhiSlices] ;
318   
319   // list of point as used in the poisson relation and the interpolation (during sum up)
320   Double_t  rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
321   for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
322     philist[k] =  gridSizePhi * k;
323     for ( Int_t i = 0 ; i < kRows ; i++ )    {
324       rlist[i] = fgkIFCRadius + i*gridSizeR ;
325       for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
326         zedlist[j]  = j * gridSizeZ ;
327       }
328     }
329   }
330
331   // ==========================================================================
332   // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
333   // Allow for different size grid spacing in R and Z directions
334   
335   const Int_t   symmetry[6] =    {1,1,-1,-1,1,1}; // shifted rod (2x), rotated clip (2x), only rod shift on OFC (1x)
336
337   // check which lookup tables are needed ---------------------------------
338
339   Bool_t needTable[6] ={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
340
341   // Shifted rods & strips
342   for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
343     if (fRodVoltShiftA[rod]!=0) needTable[0]=kTRUE;
344     if (fRodVoltShiftC[rod]!=0) needTable[0]=kTRUE;
345   }
346   for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
347     if (fRodVoltShiftA[rod]!=0) needTable[1]=kTRUE;
348     if (fRodVoltShiftC[rod]!=0) needTable[1]=kTRUE;
349   }
350   // Rotated clips
351   if (fRotatedClipVoltA[0]!=0 || fRotatedClipVoltC[0]!=0) needTable[2]=kTRUE;
352   if (fRotatedClipVoltA[1]!=0 || fRotatedClipVoltC[1]!=0) needTable[3]=kTRUE;
353  
354   // shifted Copper rods 
355   for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
356     if (fCopperRodShiftA[rod]!=0) needTable[4]=kTRUE;
357     if (fCopperRodShiftC[rod]!=0) needTable[4]=kTRUE;
358   }
359   // shifted Copper rods 
360   for ( Int_t rod = 18; rod < 36 ; rod++ ) {
361     if (fCopperRodShiftA[rod]!=0) needTable[5]=kTRUE;
362     if (fCopperRodShiftC[rod]!=0) needTable[5]=kTRUE;
363   }
364
365
366   // reserve the arrays for the basic lookup tables ----------------------
367   if (needTable[0] && !fInitLookUpBasic[0]) {
368     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
369       fLookUpBasic1ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
370       fLookUpBasic1EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
371       fLookUpBasic1DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
372       // will be deleted in destructor
373     }
374   }
375   if (needTable[1] && !fInitLookUpBasic[1]) {
376     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
377       fLookUpBasic2ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
378       fLookUpBasic2EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
379       fLookUpBasic2DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
380       // will be deleted in destructor
381     }
382   }
383   if (needTable[2] && !fInitLookUpBasic[2]) {
384     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
385       fLookUpBasic3ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
386       fLookUpBasic3EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
387       fLookUpBasic3DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
388       // will be deleted in destructor
389     }
390   }
391   if (needTable[3] && !fInitLookUpBasic[3]) {
392     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
393       fLookUpBasic4ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
394       fLookUpBasic4EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
395       fLookUpBasic4DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
396       // will be deleted in destructor
397     }
398   }
399   if (needTable[4] && !fInitLookUpBasic[4]) {
400     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
401       fLookUpBasic5ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
402       fLookUpBasic5EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
403       fLookUpBasic5DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
404       // will be deleted in destructor
405     }
406   }
407   if (needTable[5] && !fInitLookUpBasic[5]) {
408     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
409       fLookUpBasic6ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
410       fLookUpBasic6EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
411       fLookUpBasic6DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
412       // will be deleted in destructor
413     }
414   }
415  
416   // Set bondaries and solve Poisson's equation --------------------------
417  
418   for (Int_t look=0; look<6; look++) {
419    
420     Float_t inner18[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } ;  
421     Float_t outer18[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } ; 
422   
423     if (needTable[look] && !fInitLookUpBasic[look]) {
424
425       // Specify which rod is the reference; other distortions calculated by summing over multiple rotations of refrence
426       // Note: the interpolation later on depends on it!! Do not change if not really needed!
427       if (look==0) {
428         AliInfo(Form("IFC - ROD&Strip shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
429         inner18[0] = 1;  
430       } else if (look==1) {
431         AliInfo(Form("OFC - ROD&Strip shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
432         outer18[0] = 1;  
433       } else if (look==2) {
434         AliInfo(Form("IFC - Clip rot. : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
435         inner18[0] = 1; 
436       } else if (look==3) {
437         AliInfo(Form("OFC - Clip rot. : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
438         outer18[0] = 1;  
439       } else if (look==4) {
440         AliInfo(Form("IFC - CopperRod shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
441         inner18[0] = 1;  
442       } else if (look==5) {
443         AliInfo(Form("OFC - CopperRod shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
444         outer18[0] = 1;  
445       }
446       // Build potentials on boundary stripes for a rotated clip (SYMMETRY==-1) or a shifted rod (SYMMETRY==1)
447       if (look<4) {
448         // in these cases, the strips follow the actual rod shift (linear interpolation between the rods)
449         for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
450           Int_t slices = kPhiSlicesPerSector ;
451           Int_t ipoint = k/slices  ;  
452           innerList[k] = (((ipoint+1)*slices-k)*inner18[ipoint]-(k-ipoint*slices)*inner18[ipoint+1])/slices ;
453           outerList[k] = (((ipoint+1)*slices-k)*outer18[ipoint]-(k-ipoint*slices)*outer18[ipoint+1])/slices ;
454           if ( (k%slices) == 0 && symmetry[look] == -1 ) { innerList[k] = 0.0 ; outerList[k] = 0.0 ; } 
455           // above, force Zero if Anti-Sym
456         } 
457       } else if (look==4){
458         // in this case, we assume that the strips stay at the correct position, but the rods move
459         // the distortion is then much more localized around the rod (indicated by real data)
460         for ( Int_t k = 0 ; k < kPhiSlices ; k++ )
461           innerList[k] = outerList[k] = 0;
462         innerList[0] = inner18[0];     // point at rod
463         innerList[0] = inner18[0]/4*3;  // point close to rod (educated guess)
464       } else if (look==5){
465         // in this case, we assume that the strips stay at the correct position, but the copper plated OFC-rods move
466         // the distortion is then much more localized around the rod (indicated by real data)
467         for ( Int_t k = 0 ; k < kPhiSlices ; k++ )
468           innerList[k] = outerList[k] = 0;
469         outerList[0] = outer18[0];   // point at rod
470         outerList[1] = outer18[0]/4; // point close to rod (educated-`guessed` scaling)
471       }
472
473       // Fill arrays with initial conditions.  V on the boundary and Charge in the volume.
474       for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
475         TMatrixD &arrayV    =  *arrayofArrayV[k] ;
476         TMatrixD &charge    =  *arrayofCharge[k] ;
477         for ( Int_t i = 0 ; i < kRows ; i++ )    {
478           for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
479             arrayV(i,j) = 0.0 ; 
480             charge(i,j) = 0.0 ;
481             if ( i == 0 )       arrayV(i,j) = innerList[k] ; 
482             if ( i == kRows-1 ) arrayV(i,j) = outerList[k] ; 
483           }
484         }      
485         // no charge in the volume
486         for ( Int_t i = 1 ; i < kRows-1 ; i++ )  { 
487           for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
488             charge(i,j)  =  0.0 ;
489           }
490         }
491       }      
492      
493       // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
494       // Allow for different size grid spacing in R and Z directions
495       if (look==0) {
496         PoissonRelaxation3D( arrayofArrayV, arrayofCharge, 
497                              fLookUpBasic1ErOverEz, fLookUpBasic1EphiOverEz, fLookUpBasic1DeltaEz,
498                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[0]) ;
499         AliInfo("IFC - ROD&Strip shift : done ");
500       } else if (look==1) {
501         PoissonRelaxation3D( arrayofArrayV, arrayofCharge, 
502                              fLookUpBasic2ErOverEz, fLookUpBasic2EphiOverEz, fLookUpBasic2DeltaEz,
503                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[1]) ;
504         
505         AliInfo("OFC - ROD&Strip shift : done ");
506       } else if (look==2) {
507         PoissonRelaxation3D( arrayofArrayV, arrayofCharge, 
508                              fLookUpBasic3ErOverEz, fLookUpBasic3EphiOverEz, fLookUpBasic3DeltaEz,
509                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[2]) ;
510         AliInfo("IFC - Clip rot. : done ");
511       } else if (look==3) {
512         PoissonRelaxation3D( arrayofArrayV, arrayofCharge, 
513                              fLookUpBasic4ErOverEz, fLookUpBasic4EphiOverEz, fLookUpBasic4DeltaEz,
514                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[3]) ;
515         AliInfo("OFC - Clip rot. : done ");
516       } else if (look==4) {
517         PoissonRelaxation3D( arrayofArrayV, arrayofCharge, 
518                              fLookUpBasic5ErOverEz, fLookUpBasic5EphiOverEz, fLookUpBasic5DeltaEz,
519                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[4]) ;
520         AliInfo("IFC - CopperRod shift : done ");
521       } else if (look==5) {
522         PoissonRelaxation3D( arrayofArrayV, arrayofCharge, 
523                              fLookUpBasic6ErOverEz, fLookUpBasic6EphiOverEz, fLookUpBasic6DeltaEz,
524                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[5]) ;
525         AliInfo("OFC - CopperRod shift : done ");
526       }
527       
528       fInitLookUpBasic[look] = kTRUE;
529     }
530   }
531   
532
533   // =============================================================================
534   // Create the final lookup table with corresponding ROD shifts or clip rotations
535
536   Float_t erIntegralSum   = 0.0 ;
537   Float_t ephiIntegralSum = 0.0 ;
538   Float_t ezIntegralSum   = 0.0 ;
539
540   Double_t phiPrime     = 0. ;
541   Double_t erIntegral   = 0. ;
542   Double_t ephiIntegral = 0. ;
543   Double_t ezIntegral   = 0. ;
544
545
546   AliInfo("Calculation of combined Look-up Table");
547
548   // Interpolate and sum the Basic lookup tables onto the standard grid
549   Double_t  r, phi, z ;
550   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
551     phi = fgkPhiList[k] ;
552
553     TMatrixF &erOverEz   =  *fLookUpErOverEz[k]  ;
554     TMatrixF &ephiOverEz =  *fLookUpEphiOverEz[k];
555     TMatrixF &deltaEz    =  *fLookUpDeltaEz[k]   ;
556
557     for ( Int_t i = 0 ; i < kNZ ; i++ ) {
558       z = TMath::Abs(fgkZList[i]) ;  // Symmetric solution in Z that depends only on ABS(Z)
559       for ( Int_t j = 0 ; j < kNR ; j++ ) { 
560         r = fgkRList[j] ;
561         // Interpolate basicLookup tables; once for each rod, then sum the results
562         
563         erIntegralSum   = 0.0 ;
564         ephiIntegralSum = 0.0 ;
565         ezIntegralSum   = 0.0 ;
566  
567         // SHIFTED RODS =========================================================
568
569         // IFC ROD SHIFTS +++++++++++++++++++++++++++++
570         for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
571           
572           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
573           
574           if ( fRodVoltShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
575           if ( fRodVoltShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
576
577           // Rotate to a position where we have maps and prepare to sum
578           phiPrime =  phi - rod*kPhiSlicesPerSector*gridSizePhi ;  
579
580           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi    
581
582           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
583             
584             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
585                                               rlist, zedlist, philist, fLookUpBasic1ErOverEz  );
586             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
587                                               rlist, zedlist, philist, fLookUpBasic1EphiOverEz);
588             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
589                                               rlist, zedlist, philist, fLookUpBasic1DeltaEz   );
590             
591           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
592             
593             phiPrime     = TMath::TwoPi() - phiPrime ;
594             
595             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
596                                               rlist, zedlist, philist, fLookUpBasic1ErOverEz  );
597             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
598                                               rlist, zedlist, philist, fLookUpBasic1EphiOverEz);
599             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
600                                               rlist, zedlist, philist, fLookUpBasic1DeltaEz   );
601             
602             // Flip symmetry of phi integral for symmetric boundary conditions
603             if ( symmetry[0] ==  1 ) ephiIntegral *= -1  ;    
604             // Flip symmetry of r integral if anti-symmetric boundary conditions 
605             if ( symmetry[0] == -1 ) erIntegral   *= -1  ;    
606
607           }
608
609           if ( fgkZList[i] > 0 ) {
610             erIntegralSum   += fRodVoltShiftA[rod]*erIntegral   ;
611             ephiIntegralSum += fRodVoltShiftA[rod]*ephiIntegral ;
612             ezIntegralSum   += fRodVoltShiftA[rod]*ezIntegral;
613           } else if ( fgkZList[i] < 0 ) {
614             erIntegralSum   += fRodVoltShiftC[rod]*erIntegral   ;
615             ephiIntegralSum += fRodVoltShiftC[rod]*ephiIntegral ;
616             ezIntegralSum   -= fRodVoltShiftC[rod]*ezIntegral;
617           }
618         }
619
620         // OFC ROD SHIFTS +++++++++++++++++++++++++++++
621         for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
622
623           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
624           
625           if ( fRodVoltShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
626           if ( fRodVoltShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
627
628           // Rotate to a position where we have maps and prepare to sum
629           phiPrime =  phi - (rod-18)*kPhiSlicesPerSector*gridSizePhi ;  
630                   
631           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi    
632
633           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
634             
635             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
636                                               rlist, zedlist, philist, fLookUpBasic2ErOverEz  );
637             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
638                                               rlist, zedlist, philist, fLookUpBasic2EphiOverEz);
639             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
640                                               rlist, zedlist, philist, fLookUpBasic2DeltaEz   );
641             
642           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
643             
644             phiPrime     = TMath::TwoPi() - phiPrime ;
645             
646             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
647                                               rlist, zedlist, philist, fLookUpBasic2ErOverEz  );
648             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
649                                               rlist, zedlist, philist, fLookUpBasic2EphiOverEz);
650             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
651                                               rlist, zedlist, philist, fLookUpBasic2DeltaEz   );
652             
653             // Flip symmetry of phi integral for symmetric boundary conditions
654             if ( symmetry[1] ==  1 ) ephiIntegral *= -1  ;    
655             // Flip symmetry of r integral if anti-symmetric boundary conditions 
656             if ( symmetry[1] == -1 ) erIntegral   *= -1  ;    
657
658           }
659
660           if ( fgkZList[i] > 0 ) {
661             erIntegralSum   += fRodVoltShiftA[rod]*erIntegral   ;
662             ephiIntegralSum += fRodVoltShiftA[rod]*ephiIntegral ;
663             ezIntegralSum   += fRodVoltShiftA[rod]*ezIntegral;
664           } else if ( fgkZList[i] < 0 ) {
665             erIntegralSum   += fRodVoltShiftC[rod]*erIntegral   ;
666             ephiIntegralSum += fRodVoltShiftC[rod]*ephiIntegral ;
667             ezIntegralSum   -= fRodVoltShiftC[rod]*ezIntegral;
668           }
669
670         } // rod loop - shited ROD
671
672
673         // Rotated clips =========================================================
674
675         Int_t rodIFC = 11; // resistor rod positions, IFC
676         Int_t rodOFC =  3; // resistor rod positions, OFC
677         // just one rod on IFC and OFC
678
679         // IFC ROTATED CLIP +++++++++++++++++++++++++++++
680         for ( Int_t rod = rodIFC ; rod < rodIFC+1 ; rod++ ) { // loop over 1 to keep the "ignore"-functionality
681
682           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
683           if ( fRotatedClipVoltA[0] == 0 && fgkZList[i] > 0) continue ;
684           if ( fRotatedClipVoltC[0] == 0 && fgkZList[i] < 0) continue ;
685
686           // Rotate to a position where we have maps and prepare to sum
687           phiPrime =  phi - rod*kPhiSlicesPerSector*gridSizePhi ;  
688           
689           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi    
690           
691           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
692             
693             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
694                                               rlist, zedlist, philist, fLookUpBasic3ErOverEz  );
695             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
696                                               rlist, zedlist, philist, fLookUpBasic3EphiOverEz);
697             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
698                                               rlist, zedlist, philist, fLookUpBasic3DeltaEz   );
699             
700           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
701             
702             phiPrime     = TMath::TwoPi() - phiPrime ;
703             
704             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
705                                               rlist, zedlist, philist, fLookUpBasic3ErOverEz  );
706             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
707                                               rlist, zedlist, philist, fLookUpBasic3EphiOverEz);
708             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
709                                               rlist, zedlist, philist, fLookUpBasic3DeltaEz   );
710             
711             // Flip symmetry of phi integral for symmetric boundary conditions
712             if ( symmetry[2] ==  1 ) ephiIntegral *= -1  ;    
713             // Flip symmetry of r integral if anti-symmetric boundary conditions 
714             if ( symmetry[2] == -1 ) erIntegral   *= -1  ;    
715             
716           }
717           
718           if ( fgkZList[i] > 0 ) {
719             erIntegralSum   += fRotatedClipVoltA[0]*erIntegral   ;
720             ephiIntegralSum += fRotatedClipVoltA[0]*ephiIntegral ;
721             ezIntegralSum   += fRotatedClipVoltA[0]*ezIntegral;
722           } else if ( fgkZList[i] < 0 ) {
723             erIntegralSum   += fRotatedClipVoltC[0]*erIntegral   ;
724             ephiIntegralSum += fRotatedClipVoltC[0]*ephiIntegral ;
725             ezIntegralSum   -= fRotatedClipVoltC[0]*ezIntegral;
726           }
727         }
728
729         // OFC: ROTATED CLIP  +++++++++++++++++++++++++++++
730         for ( Int_t rod = rodOFC ; rod < rodOFC+1 ; rod++ ) { // loop over 1 to keep the "ignore"-functionality
731           
732           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
733           
734           if ( fRotatedClipVoltA[1] == 0 && fgkZList[i] > 0) continue ;
735           if ( fRotatedClipVoltC[1] == 0 && fgkZList[i] < 0) continue ;
736
737           // Rotate to a position where we have maps and prepare to sum
738           phiPrime =  phi - rod*kPhiSlicesPerSector*gridSizePhi ;  
739           
740           
741           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi    
742           
743           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
744             
745             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
746                                               rlist, zedlist, philist, fLookUpBasic4ErOverEz  );
747             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
748                                               rlist, zedlist, philist, fLookUpBasic4EphiOverEz);
749             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
750                                               rlist, zedlist, philist, fLookUpBasic4DeltaEz   );
751             
752           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
753             
754             phiPrime     = TMath::TwoPi() - phiPrime ;
755             
756             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
757                                               rlist, zedlist, philist, fLookUpBasic4ErOverEz  );
758             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
759                                               rlist, zedlist, philist, fLookUpBasic4EphiOverEz);
760             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
761                                               rlist, zedlist, philist, fLookUpBasic4DeltaEz   );
762             
763             // Flip symmetry of phi integral for symmetric boundary conditions
764             if ( symmetry[3] ==  1 ) ephiIntegral *= -1  ;    
765             // Flip symmetry of r integral if anti-symmetric boundary conditions 
766             if ( symmetry[3] == -1 ) erIntegral   *= -1  ;    
767             
768           }
769           
770           if ( fgkZList[i] > 0 ) {
771             erIntegralSum   += fRotatedClipVoltA[1]*erIntegral   ;
772             ephiIntegralSum += fRotatedClipVoltA[1]*ephiIntegral ;
773             ezIntegralSum   += fRotatedClipVoltA[1]*ezIntegral;
774           } else if ( fgkZList[i] < 0 ) {
775             erIntegralSum   += fRotatedClipVoltC[1]*erIntegral   ;
776             ephiIntegralSum += fRotatedClipVoltC[1]*ephiIntegral ;
777             ezIntegralSum   -= fRotatedClipVoltC[1]*ezIntegral;
778           }
779         }
780
781         // IFC Cooper rod shift  +++++++++++++++++++++++++++++
782         for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
783           
784           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
785           
786           if ( fCopperRodShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
787           if ( fCopperRodShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
788
789           // Rotate to a position where we have maps and prepare to sum
790           phiPrime =  phi - rod*kPhiSlicesPerSector*gridSizePhi ;  
791
792           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi    
793
794           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
795             
796             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
797                                               rlist, zedlist, philist, fLookUpBasic5ErOverEz  );
798             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
799                                               rlist, zedlist, philist, fLookUpBasic5EphiOverEz);
800             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
801                                               rlist, zedlist, philist, fLookUpBasic5DeltaEz   );
802             
803           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
804             
805             phiPrime     = TMath::TwoPi() - phiPrime ;
806             
807             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
808                                               rlist, zedlist, philist, fLookUpBasic5ErOverEz  );
809             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
810                                               rlist, zedlist, philist, fLookUpBasic5EphiOverEz);
811             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
812                                               rlist, zedlist, philist, fLookUpBasic5DeltaEz   );
813             
814             // Flip symmetry of phi integral for symmetric boundary conditions
815             if ( symmetry[4] ==  1 ) ephiIntegral *= -1  ;    
816             // Flip symmetry of r integral if anti-symmetric boundary conditions 
817             if ( symmetry[4] == -1 ) erIntegral   *= -1  ;    
818
819           }
820
821           if ( fgkZList[i] > 0 ) {
822             erIntegralSum   += fCopperRodShiftA[rod]*erIntegral   ;
823             ephiIntegralSum += fCopperRodShiftA[rod]*ephiIntegral ;
824             ezIntegralSum   += fCopperRodShiftA[rod]*ezIntegral;
825           } else if ( fgkZList[i] < 0 ) {
826             erIntegralSum   += fCopperRodShiftC[rod]*erIntegral   ;
827             ephiIntegralSum += fCopperRodShiftC[rod]*ephiIntegral ;
828             ezIntegralSum   -= fCopperRodShiftC[rod]*ezIntegral;
829           }
830         }
831
832         // OFC Cooper rod shift  +++++++++++++++++++++++++++++
833         for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
834           
835           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
836           
837           if ( fCopperRodShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
838           if ( fCopperRodShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
839
840           // Rotate to a position where we have maps and prepare to sum
841           phiPrime =  phi - (rod-18)*kPhiSlicesPerSector*gridSizePhi ;  
842
843           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi    
844
845           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
846             
847             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
848                                               rlist, zedlist, philist, fLookUpBasic6ErOverEz  );
849             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
850                                               rlist, zedlist, philist, fLookUpBasic6EphiOverEz);
851             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
852                                               rlist, zedlist, philist, fLookUpBasic6DeltaEz   );
853             
854           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
855             
856             phiPrime     = TMath::TwoPi() - phiPrime ;
857             
858             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices, 
859                                               rlist, zedlist, philist, fLookUpBasic6ErOverEz  );
860             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
861                                               rlist, zedlist, philist, fLookUpBasic6EphiOverEz);
862             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
863                                               rlist, zedlist, philist, fLookUpBasic6DeltaEz   );
864             
865             // Flip symmetry of phi integral for symmetric boundary conditions
866             if ( symmetry[5] ==  1 ) ephiIntegral *= -1  ;    
867             // Flip symmetry of r integral if anti-symmetric boundary conditions 
868             if ( symmetry[5] == -1 ) erIntegral   *= -1  ;    
869
870           }
871
872           if ( fgkZList[i] > 0 ) {
873             erIntegralSum   += fCopperRodShiftA[rod]*erIntegral   ;
874             ephiIntegralSum += fCopperRodShiftA[rod]*ephiIntegral ;
875             ezIntegralSum   += fCopperRodShiftA[rod]*ezIntegral;
876           } else if ( fgkZList[i] < 0 ) {
877             erIntegralSum   += fCopperRodShiftC[rod]*erIntegral   ;
878             ephiIntegralSum += fCopperRodShiftC[rod]*ephiIntegral ;
879             ezIntegralSum   -= fCopperRodShiftC[rod]*ezIntegral;
880           }
881         }
882
883         // put the sum into the final lookup table
884         erOverEz(j,i)   =  erIntegralSum;
885         ephiOverEz(j,i) =  ephiIntegralSum;
886         deltaEz(j,i)    =  ezIntegralSum;
887         
888         //      if (j==1) printf("%lf %lf %lf: %lf %lf %lf\n",r, phi, z, erIntegralSum,ephiIntegralSum,ezIntegralSum);
889  
890       } // endo r loop
891     } // end of z loop
892   } // end of phi loop
893
894
895   // clear the temporary arrays lists
896   for ( Int_t k = 0 ; k < kPhiSlices ; k++ )  {
897     delete arrayofArrayV[k];
898     delete arrayofCharge[k];
899   }
900  
901   AliInfo(" done");
902   fInitLookUp = kTRUE;
903
904 }
905
906 void AliTPCFCVoltError3D::Print(const Option_t* option) const {
907   //
908   // Print function to check the settings of the Rod shifts and the rotated clips
909   // option=="a" prints the C0 and C1 coefficents for calibration purposes
910   //
911
912   TString opt = option; opt.ToLower();
913   printf("%s\n",GetTitle());
914   printf(" - ROD shifts  (residual voltage settings): 40V correspond to 1mm of shift\n");
915   Int_t count = 0;
916   printf("  : A-side - (Rod & Strips) \n     "); 
917   for (Int_t i=0; i<36;i++) {
918     if (fRodVoltShiftA[i]!=0) {
919       printf("Rod%2d:%3.1fV ",i,fRodVoltShiftA[i]);
920       count++;
921     }
922     if (count==0&&i==35) 
923       printf("-> all at 0 V\n");
924     else if (i==35){
925       printf("\n");
926       count=0;
927     }
928   } 
929   printf("  : C-side - (Rod & Strips) \n     "); 
930   for (Int_t i=0; i<36;i++) {
931     if (fRodVoltShiftC[i]!=0) {
932       printf("Rod%2d:%3.1fV ",i,fRodVoltShiftC[i]);
933       count++;
934     }
935     if (count==0&&i==35) 
936       printf("-> all at 0 V\n");
937     else if (i==35){
938       printf("\n");
939       count=0;
940     }
941   } 
942  
943   printf(" - Rotated clips (residual voltage settings): 40V correspond to 1mm of 'offset'\n");
944   if (fRotatedClipVoltA[0]!=0) {   printf("     A side (IFC): HV Rod: %3.1f V \n",fRotatedClipVoltA[0]); count++; }
945   if (fRotatedClipVoltA[1]!=0) {   printf("     A side (OFC): HV Rod: %3.1f V \n",fRotatedClipVoltA[1]); count++; }
946   if (fRotatedClipVoltC[0]!=0) {   printf("     C side (IFC): HV Rod: %3.1f V \n",fRotatedClipVoltC[0]); count++; }
947   if (fRotatedClipVoltC[1]!=0) {   printf("     C side (OFC): HV Rod: %3.1f V \n",fRotatedClipVoltC[1]); count++; }
948   if (count==0) 
949     printf("     -> no rotated clips \n");
950
951   count=0;
952   printf(" - Copper ROD shifts (without strips):\n");
953   printf("  : A-side - (Copper Rod shift) \n     "); 
954   for (Int_t i=0; i<36;i++) {
955     if (fCopperRodShiftA[i]!=0) {
956       printf("Rod%2d:%3.1fV ",i,fCopperRodShiftA[i]);
957       count++;
958     }
959     if (count==0&&i==35) 
960       printf("-> all at 0 V\n");
961     else if (i==35){
962       printf("\n");
963       count=0;
964     }
965   } 
966   printf("  : C-side - (Copper Rod shift) \n     "); 
967   for (Int_t i=0; i<36;i++) {
968     if (fCopperRodShiftC[i]!=0) {
969       printf("Rod%2d:%3.1fV ",i,fCopperRodShiftC[i]);
970       count++;
971     }
972     if (count==0&&i==35) 
973       printf("-> all at 0 V\n");
974     else if (i==35){
975       printf("\n");
976       count=0;
977     }
978   } 
979
980   if (opt.Contains("a")) { // Print all details
981     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
982     printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
983   }
984
985   if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitFCVoltError3D() ...");
986
987 }