TPC/Base/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("cAliTPCVoltError3D","cAliTPCVoltError3D",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 }