TPC/TPCbase/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
 192 Bool_t AliTPCFCVoltError3D::AddCorrectionCompact(AliTPCCorrection* corr, Double_t weight){
 193   //
 194   // Add correction  and make them compact
 195   // Assumptions:
 196   //  - origin of distortion/correction are additive
 197   //  - only correction ot the same type supported ()
 198   if (corr==NULL) {
 199     AliError("Zerro pointer - correction");
 200     return kFALSE;
 201   }
 202   AliTPCFCVoltError3D * corrC = dynamic_cast<AliTPCFCVoltError3D *>(corr);
 203   if (corrC == NULL)  {
 204     AliError(TString::Format("Inconsistent class types: %s\%s",IsA()->GetName(),corr->IsA()->GetName()).Data());
 205     return kFALSE;
 206   }
 207   //
 208   for (Int_t isec=0; isec<36; isec++){
 209     fRodVoltShiftA[isec]+= weight*corrC->fRodVoltShiftA[isec];      // Rod (&strips) shift in Volt (40V~1mm)
 210     fRodVoltShiftC[isec]+=weight*corrC->fRodVoltShiftC[isec];      // Rod (&strips) shift in Volt (40V~1mm)
 211     fCopperRodShiftA[isec]+=weight*corrC->fCopperRodShiftA[isec];    // only Rod shift
 212     fCopperRodShiftC[isec]+=weight*corrC->fCopperRodShiftC[isec];    // only Rod shift
 213   }
 214   for (Int_t isec=0; isec<2; isec++){
 215     fRotatedClipVoltA[isec]+= weight*corrC->fRotatedClipVoltA[isec];    // rotated clips at HV rod
 216     fRotatedClipVoltC[isec]+= weight*corrC-> fRotatedClipVoltC[isec];    // rotated clips at HV rod
 217   }
 218
 219   return kTRUE;
 220 }
 221
 222
 223
 224 void AliTPCFCVoltError3D::Init() {
 225   //
 226   // Initialization funtion
 227   //
 228
 229   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 230   if (!magF) AliError("Magneticd field - not initialized");
 231   Double_t bzField = magF->SolenoidField()/10.; //field in T
 232   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 233   if (!param) AliError("Parameters - not initialized");
 234   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
 235   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 236   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 237   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 238   SetOmegaTauT1T2(wt,fT1,fT2);
 239
 240   if (!fInitLookUp) InitFCVoltError3D();
 241 }
 242
 243 void AliTPCFCVoltError3D::Update(const TTimeStamp &/*timeStamp*/) {
 244   //
 245   // Update function
 246   //
 247   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 248   if (!magF) AliError("Magneticd field - not initialized");
 249   Double_t bzField = magF->SolenoidField()/10.; //field in T
 250   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 251   if (!param) AliError("Parameters - not initialized");
 252   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]  // From dataBase: to be updated: per second (ideally)
 253   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 254   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 255   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 256   SetOmegaTauT1T2(wt,fT1,fT2);
 257
 258
 259 }
 260
 261
 262
 263 void AliTPCFCVoltError3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
 264   //
 265   // Calculates the correction due e.g. residual voltage errors on the TPC boundaries
 266   //
 267   const Double_t kEpsilon=Double_t(FLT_MIN);
 268
 269   if (!fInitLookUp) {
 270     AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
 271     InitFCVoltError3D();
 272   }
 273
 274   static Bool_t forceInit=kTRUE; //temporary needed for back compatibility with old OCDB entries
 275   if (forceInit &&fLookUpErOverEz[0]){
 276     if (TMath::Abs(fLookUpErOverEz[0]->Sum())<kEpsilon){//temporary needed for back compatibility with old OCDB entries
 277       ForceInitFCVoltError3D();
 278     }
 279     forceInit=kFALSE;
 280   }
 281
 282
 283   Int_t   order     = 1 ;               // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
 284                                         // note that the poisson solution was linearly mirroed on this grid!
 285   Float_t intEr, intEphi, intDeltaEz;
 286   Float_t r, phi, z ;
 287   Int_t    sign;
 288
 289   r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
 290   phi    =  TMath::ATan2(x[1],x[0]) ;
 291   if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
 292   z      =  x[2] ;                                         // Create temporary copy of x[2]
 293
 294   if ( (roc%36) < 18 ) {
 295     sign =  1;       // (TPC A side)
 296   } else {
 297     sign = -1;       // (TPC C side)
 298   }
 299
 300   if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
 301   if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
 302
 303
 304   if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
 305     AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
 306
 307   // Get the Er and Ephi field integrals plus the integral over DeltaEz
 308   intEr      = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 309                                   fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz  );
 310   intEphi    = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 311                                   fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz  );
 312   intDeltaEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 313                                   fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz  );
 314
 315   //  printf("%lf %lf %lf\n",intEr,intEphi,intDeltaEz);
 316
 317   // Calculate distorted position
 318   if ( r > 0.0 ) {
 319     phi =  phi + ( fC0*intEphi - fC1*intEr ) / r;
 320     r   =  r   + ( fC0*intEr   + fC1*intEphi );
 321   }
 322
 323   // Calculate correction in cartesian coordinates
 324   dx[0] = r * TMath::Cos(phi) - x[0];
 325   dx[1] = r * TMath::Sin(phi) - x[1];
 326   dx[2] = intDeltaEz;  // z distortion - (internally scaled with driftvelocity dependency
 327                        // on the Ez field plus the actual ROC misalignment (if set TRUE)
 328
 329 }
 330
 331 void AliTPCFCVoltError3D::InitFCVoltError3D() {
 332   //
 333   // Initialization of the Lookup table which contains the solutions of the
 334   // Dirichlet boundary problem
 335   // Calculation of the single 3D-Poisson solver is done just if needed
 336   // (see basic lookup tables in header file)
 337   //
 338
 339   const Int_t   order       =    1  ;  // Linear interpolation = 1, Quadratic = 2
 340   const Float_t gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
 341   const Float_t gridSizeZ   =  fgkTPCZ0 / (kColumns-1) ;
 342   const Float_t gridSizePhi =  TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);
 343
 344   // temporary arrays to create the boundary conditions
 345   TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ;
 346   for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
 347     arrayofArrayV[k]     =   new TMatrixD(kRows,kColumns) ;
 348     arrayofCharge[k]     =   new TMatrixD(kRows,kColumns) ;
 349   }
 350   Double_t  innerList[kPhiSlices], outerList[kPhiSlices] ;
 351
 352   // list of point as used in the poisson relation and the interpolation (during sum up)
 353   Double_t  rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
 354   for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
 355     philist[k] =  gridSizePhi * k;
 356     for ( Int_t i = 0 ; i < kRows ; i++ )    {
 357       rlist[i] = fgkIFCRadius + i*gridSizeR ;
 358       for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
 359         zedlist[j]  = j * gridSizeZ ;
 360       }
 361     }
 362   }
 363
 364   // ==========================================================================
 365   // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
 366   // Allow for different size grid spacing in R and Z directions
 367
 368   const Int_t   symmetry[6] =    {1,1,-1,-1,1,1}; // shifted rod (2x), rotated clip (2x), only rod shift on OFC (1x)
 369
 370   // check which lookup tables are needed ---------------------------------
 371
 372   Bool_t needTable[6] ={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
 373
 374   // Shifted rods & strips
 375   for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
 376     if (fRodVoltShiftA[rod]!=0) needTable[0]=kTRUE;
 377     if (fRodVoltShiftC[rod]!=0) needTable[0]=kTRUE;
 378   }
 379   for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
 380     if (fRodVoltShiftA[rod]!=0) needTable[1]=kTRUE;
 381     if (fRodVoltShiftC[rod]!=0) needTable[1]=kTRUE;
 382   }
 383   // Rotated clips
 384   if (fRotatedClipVoltA[0]!=0 || fRotatedClipVoltC[0]!=0) needTable[2]=kTRUE;
 385   if (fRotatedClipVoltA[1]!=0 || fRotatedClipVoltC[1]!=0) needTable[3]=kTRUE;
 386
 387   // shifted Copper rods
 388   for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
 389     if (fCopperRodShiftA[rod]!=0) needTable[4]=kTRUE;
 390     if (fCopperRodShiftC[rod]!=0) needTable[4]=kTRUE;
 391   }
 392   // shifted Copper rods
 393   for ( Int_t rod = 18; rod < 36 ; rod++ ) {
 394     if (fCopperRodShiftA[rod]!=0) needTable[5]=kTRUE;
 395     if (fCopperRodShiftC[rod]!=0) needTable[5]=kTRUE;
 396   }
 397
 398
 399   // reserve the arrays for the basic lookup tables ----------------------
 400   if (needTable[0] && !fInitLookUpBasic[0]) {
 401     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
 402       fLookUpBasic1ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
 403       fLookUpBasic1EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
 404       fLookUpBasic1DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
 405       // will be deleted in destructor
 406     }
 407   }
 408   if (needTable[1] && !fInitLookUpBasic[1]) {
 409     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
 410       fLookUpBasic2ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
 411       fLookUpBasic2EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
 412       fLookUpBasic2DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
 413       // will be deleted in destructor
 414     }
 415   }
 416   if (needTable[2] && !fInitLookUpBasic[2]) {
 417     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
 418       fLookUpBasic3ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
 419       fLookUpBasic3EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
 420       fLookUpBasic3DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
 421       // will be deleted in destructor
 422     }
 423   }
 424   if (needTable[3] && !fInitLookUpBasic[3]) {
 425     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
 426       fLookUpBasic4ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
 427       fLookUpBasic4EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
 428       fLookUpBasic4DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
 429       // will be deleted in destructor
 430     }
 431   }
 432   if (needTable[4] && !fInitLookUpBasic[4]) {
 433     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
 434       fLookUpBasic5ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
 435       fLookUpBasic5EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
 436       fLookUpBasic5DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
 437       // will be deleted in destructor
 438     }
 439   }
 440   if (needTable[5] && !fInitLookUpBasic[5]) {
 441     for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {   // Possibly make an extra table to be used for 0 == 360
 442       fLookUpBasic6ErOverEz[k]   =   new TMatrixD(kRows,kColumns);
 443       fLookUpBasic6EphiOverEz[k] =   new TMatrixD(kRows,kColumns);
 444       fLookUpBasic6DeltaEz[k]    =   new TMatrixD(kRows,kColumns);
 445       // will be deleted in destructor
 446     }
 447   }
 448
 449   // Set bondaries and solve Poisson's equation --------------------------
 450
 451   for (Int_t look=0; look<6; look++) {
 452
 453     Float_t inner18[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } ;
 454     Float_t outer18[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } ;
 455
 456     if (needTable[look] && !fInitLookUpBasic[look]) {
 457
 458       // Specify which rod is the reference; other distortions calculated by summing over multiple rotations of refrence
 459       // Note: the interpolation later on depends on it!! Do not change if not really needed!
 460       if (look==0) {
 461         AliInfo(Form("IFC - ROD&Strip shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
 462         inner18[0] = 1;
 463       } else if (look==1) {
 464         AliInfo(Form("OFC - ROD&Strip shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
 465         outer18[0] = 1;
 466       } else if (look==2) {
 467         AliInfo(Form("IFC - Clip rot. : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
 468         inner18[0] = 1;
 469       } else if (look==3) {
 470         AliInfo(Form("OFC - Clip rot. : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
 471         outer18[0] = 1;
 472       } else if (look==4) {
 473         AliInfo(Form("IFC - CopperRod shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
 474         inner18[0] = 1;
 475       } else if (look==5) {
 476         AliInfo(Form("OFC - CopperRod shift : Filling table (~ %d sec)",3*(int)(kPhiSlices/5)));
 477         outer18[0] = 1;
 478       }
 479       // Build potentials on boundary stripes for a rotated clip (SYMMETRY==-1) or a shifted rod (SYMMETRY==1)
 480       if (look<4) {
 481         // in these cases, the strips follow the actual rod shift (linear interpolation between the rods)
 482         for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
 483           Int_t slices = kPhiSlicesPerSector ;
 484           Int_t ipoint = k/slices  ;
 485           innerList[k] = (((ipoint+1)*slices-k)*inner18[ipoint]-(k-ipoint*slices)*inner18[ipoint+1])/slices ;
 486           outerList[k] = (((ipoint+1)*slices-k)*outer18[ipoint]-(k-ipoint*slices)*outer18[ipoint+1])/slices ;
 487           if ( (k%slices) == 0 && symmetry[look] == -1 ) { innerList[k] = 0.0 ; outerList[k] = 0.0 ; }
 488           // above, force Zero if Anti-Sym
 489         }
 490       } else if (look==4){
 491         // in this case, we assume that the strips stay at the correct position, but the rods move
 492         // the distortion is then much more localized around the rod (indicated by real data)
 493         for ( Int_t k = 0 ; k < kPhiSlices ; k++ )
 494           innerList[k] = outerList[k] = 0;
 495         innerList[0] = inner18[0];     // point at rod
 496         innerList[0] = inner18[0]/4*3;  // point close to rod (educated guess)
 497       } else if (look==5){
 498         // in this case, we assume that the strips stay at the correct position, but the copper plated OFC-rods move
 499         // the distortion is then much more localized around the rod (indicated by real data)
 500         for ( Int_t k = 0 ; k < kPhiSlices ; k++ )
 501           innerList[k] = outerList[k] = 0;
 502         outerList[0] = outer18[0];   // point at rod
 503         outerList[1] = outer18[0]/4; // point close to rod (educated-`guessed` scaling)
 504       }
 505
 506       // Fill arrays with initial conditions.  V on the boundary and Charge in the volume.
 507       for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
 508         TMatrixD &arrayV    =  *arrayofArrayV[k] ;
 509         TMatrixD &charge    =  *arrayofCharge[k] ;
 510         for ( Int_t i = 0 ; i < kRows ; i++ )    {
 511           for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
 512             arrayV(i,j) = 0.0 ;
 513             charge(i,j) = 0.0 ;
 514             if ( i == 0 )       arrayV(i,j) = innerList[k] ;
 515             if ( i == kRows-1 ) arrayV(i,j) = outerList[k] ;
 516           }
 517         }
 518         // no charge in the volume
 519         for ( Int_t i = 1 ; i < kRows-1 ; i++ )  {
 520           for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
 521             charge(i,j)  =  0.0 ;
 522           }
 523         }
 524       }
 525
 526       // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
 527       // Allow for different size grid spacing in R and Z directions
 528       if (look==0) {
 529         PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
 530                              fLookUpBasic1ErOverEz, fLookUpBasic1EphiOverEz, fLookUpBasic1DeltaEz,
 531                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[0]) ;
 532         AliInfo("IFC - ROD&Strip shift : done ");
 533       } else if (look==1) {
 534         PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
 535                              fLookUpBasic2ErOverEz, fLookUpBasic2EphiOverEz, fLookUpBasic2DeltaEz,
 536                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[1]) ;
 537
 538         AliInfo("OFC - ROD&Strip shift : done ");
 539       } else if (look==2) {
 540         PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
 541                              fLookUpBasic3ErOverEz, fLookUpBasic3EphiOverEz, fLookUpBasic3DeltaEz,
 542                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[2]) ;
 543         AliInfo("IFC - Clip rot. : done ");
 544       } else if (look==3) {
 545         PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
 546                              fLookUpBasic4ErOverEz, fLookUpBasic4EphiOverEz, fLookUpBasic4DeltaEz,
 547                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[3]) ;
 548         AliInfo("OFC - Clip rot. : done ");
 549       } else if (look==4) {
 550         PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
 551                              fLookUpBasic5ErOverEz, fLookUpBasic5EphiOverEz, fLookUpBasic5DeltaEz,
 552                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[4]) ;
 553         AliInfo("IFC - CopperRod shift : done ");
 554       } else if (look==5) {
 555         PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
 556                              fLookUpBasic6ErOverEz, fLookUpBasic6EphiOverEz, fLookUpBasic6DeltaEz,
 557                              kRows, kColumns, kPhiSlices, gridSizePhi, kIterations, symmetry[5]) ;
 558         AliInfo("OFC - CopperRod shift : done ");
 559       }
 560
 561       fInitLookUpBasic[look] = kTRUE;
 562     }
 563   }
 564
 565
 566   // =============================================================================
 567   // Create the final lookup table with corresponding ROD shifts or clip rotations
 568
 569   Float_t erIntegralSum   = 0.0 ;
 570   Float_t ephiIntegralSum = 0.0 ;
 571   Float_t ezIntegralSum   = 0.0 ;
 572
 573   Double_t phiPrime     = 0. ;
 574   Double_t erIntegral   = 0. ;
 575   Double_t ephiIntegral = 0. ;
 576   Double_t ezIntegral   = 0. ;
 577
 578
 579   AliInfo("Calculation of combined Look-up Table");
 580
 581   // Interpolate and sum the Basic lookup tables onto the standard grid
 582   Double_t  r, phi, z ;
 583   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 584     phi = fgkPhiList[k] ;
 585
 586     TMatrixF &erOverEz   =  *fLookUpErOverEz[k]  ;
 587     TMatrixF &ephiOverEz =  *fLookUpEphiOverEz[k];
 588     TMatrixF &deltaEz    =  *fLookUpDeltaEz[k]   ;
 589
 590     for ( Int_t i = 0 ; i < kNZ ; i++ ) {
 591       z = TMath::Abs(fgkZList[i]) ;  // Symmetric solution in Z that depends only on ABS(Z)
 592       for ( Int_t j = 0 ; j < kNR ; j++ ) {
 593         r = fgkRList[j] ;
 594         // Interpolate basicLookup tables; once for each rod, then sum the results
 595
 596         erIntegralSum   = 0.0 ;
 597         ephiIntegralSum = 0.0 ;
 598         ezIntegralSum   = 0.0 ;
 599
 600         // SHIFTED RODS =========================================================
 601
 602         // IFC ROD SHIFTS +++++++++++++++++++++++++++++
 603         for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
 604
 605           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
 606
 607           if ( fRodVoltShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
 608           if ( fRodVoltShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
 609
 610           // Rotate to a position where we have maps and prepare to sum
 611           phiPrime =  phi - rod*kPhiSlicesPerSector*gridSizePhi ;
 612
 613           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi
 614
 615           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
 616
 617             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 618                                               rlist, zedlist, philist, fLookUpBasic1ErOverEz  );
 619             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 620                                               rlist, zedlist, philist, fLookUpBasic1EphiOverEz);
 621             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 622                                               rlist, zedlist, philist, fLookUpBasic1DeltaEz   );
 623
 624           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
 625
 626             phiPrime     = TMath::TwoPi() - phiPrime ;
 627
 628             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 629                                               rlist, zedlist, philist, fLookUpBasic1ErOverEz  );
 630             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 631                                               rlist, zedlist, philist, fLookUpBasic1EphiOverEz);
 632             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 633                                               rlist, zedlist, philist, fLookUpBasic1DeltaEz   );
 634
 635             // Flip symmetry of phi integral for symmetric boundary conditions
 636             if ( symmetry[0] ==  1 ) ephiIntegral *= -1  ;
 637             // Flip symmetry of r integral if anti-symmetric boundary conditions
 638             if ( symmetry[0] == -1 ) erIntegral   *= -1  ;
 639
 640           }
 641
 642           if ( fgkZList[i] > 0 ) {
 643             erIntegralSum   += fRodVoltShiftA[rod]*erIntegral   ;
 644             ephiIntegralSum += fRodVoltShiftA[rod]*ephiIntegral ;
 645             ezIntegralSum   += fRodVoltShiftA[rod]*ezIntegral;
 646           } else if ( fgkZList[i] < 0 ) {
 647             erIntegralSum   += fRodVoltShiftC[rod]*erIntegral   ;
 648             ephiIntegralSum += fRodVoltShiftC[rod]*ephiIntegral ;
 649             ezIntegralSum   -= fRodVoltShiftC[rod]*ezIntegral;
 650           }
 651         }
 652
 653         // OFC ROD SHIFTS +++++++++++++++++++++++++++++
 654         for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
 655
 656           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
 657
 658           if ( fRodVoltShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
 659           if ( fRodVoltShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
 660
 661           // Rotate to a position where we have maps and prepare to sum
 662           phiPrime =  phi - (rod-18)*kPhiSlicesPerSector*gridSizePhi ;
 663
 664           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi
 665
 666           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
 667
 668             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 669                                               rlist, zedlist, philist, fLookUpBasic2ErOverEz  );
 670             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 671                                               rlist, zedlist, philist, fLookUpBasic2EphiOverEz);
 672             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 673                                               rlist, zedlist, philist, fLookUpBasic2DeltaEz   );
 674
 675           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
 676
 677             phiPrime     = TMath::TwoPi() - phiPrime ;
 678
 679             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 680                                               rlist, zedlist, philist, fLookUpBasic2ErOverEz  );
 681             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 682                                               rlist, zedlist, philist, fLookUpBasic2EphiOverEz);
 683             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 684                                               rlist, zedlist, philist, fLookUpBasic2DeltaEz   );
 685
 686             // Flip symmetry of phi integral for symmetric boundary conditions
 687             if ( symmetry[1] ==  1 ) ephiIntegral *= -1  ;
 688             // Flip symmetry of r integral if anti-symmetric boundary conditions
 689             if ( symmetry[1] == -1 ) erIntegral   *= -1  ;
 690
 691           }
 692
 693           if ( fgkZList[i] > 0 ) {
 694             erIntegralSum   += fRodVoltShiftA[rod]*erIntegral   ;
 695             ephiIntegralSum += fRodVoltShiftA[rod]*ephiIntegral ;
 696             ezIntegralSum   += fRodVoltShiftA[rod]*ezIntegral;
 697           } else if ( fgkZList[i] < 0 ) {
 698             erIntegralSum   += fRodVoltShiftC[rod]*erIntegral   ;
 699             ephiIntegralSum += fRodVoltShiftC[rod]*ephiIntegral ;
 700             ezIntegralSum   -= fRodVoltShiftC[rod]*ezIntegral;
 701           }
 702
 703         } // rod loop - shited ROD
 704
 705
 706         // Rotated clips =========================================================
 707
 708         Int_t rodIFC = 11; // resistor rod positions, IFC
 709         Int_t rodOFC =  3; // resistor rod positions, OFC
 710         // just one rod on IFC and OFC
 711
 712         // IFC ROTATED CLIP +++++++++++++++++++++++++++++
 713         for ( Int_t rod = rodIFC ; rod < rodIFC+1 ; rod++ ) { // loop over 1 to keep the "ignore"-functionality
 714
 715           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
 716           if ( fRotatedClipVoltA[0] == 0 && fgkZList[i] > 0) continue ;
 717           if ( fRotatedClipVoltC[0] == 0 && fgkZList[i] < 0) continue ;
 718
 719           // Rotate to a position where we have maps and prepare to sum
 720           phiPrime =  phi - rod*kPhiSlicesPerSector*gridSizePhi ;
 721
 722           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi
 723
 724           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
 725
 726             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 727                                               rlist, zedlist, philist, fLookUpBasic3ErOverEz  );
 728             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 729                                               rlist, zedlist, philist, fLookUpBasic3EphiOverEz);
 730             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 731                                               rlist, zedlist, philist, fLookUpBasic3DeltaEz   );
 732
 733           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
 734
 735             phiPrime     = TMath::TwoPi() - phiPrime ;
 736
 737             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 738                                               rlist, zedlist, philist, fLookUpBasic3ErOverEz  );
 739             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 740                                               rlist, zedlist, philist, fLookUpBasic3EphiOverEz);
 741             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 742                                               rlist, zedlist, philist, fLookUpBasic3DeltaEz   );
 743
 744             // Flip symmetry of phi integral for symmetric boundary conditions
 745             if ( symmetry[2] ==  1 ) ephiIntegral *= -1  ;
 746             // Flip symmetry of r integral if anti-symmetric boundary conditions
 747             if ( symmetry[2] == -1 ) erIntegral   *= -1  ;
 748
 749           }
 750
 751           if ( fgkZList[i] > 0 ) {
 752             erIntegralSum   += fRotatedClipVoltA[0]*erIntegral   ;
 753             ephiIntegralSum += fRotatedClipVoltA[0]*ephiIntegral ;
 754             ezIntegralSum   += fRotatedClipVoltA[0]*ezIntegral;
 755           } else if ( fgkZList[i] < 0 ) {
 756             erIntegralSum   += fRotatedClipVoltC[0]*erIntegral   ;
 757             ephiIntegralSum += fRotatedClipVoltC[0]*ephiIntegral ;
 758             ezIntegralSum   -= fRotatedClipVoltC[0]*ezIntegral;
 759           }
 760         }
 761
 762         // OFC: ROTATED CLIP  +++++++++++++++++++++++++++++
 763         for ( Int_t rod = rodOFC ; rod < rodOFC+1 ; rod++ ) { // loop over 1 to keep the "ignore"-functionality
 764
 765           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
 766
 767           if ( fRotatedClipVoltA[1] == 0 && fgkZList[i] > 0) continue ;
 768           if ( fRotatedClipVoltC[1] == 0 && fgkZList[i] < 0) continue ;
 769
 770           // Rotate to a position where we have maps and prepare to sum
 771           phiPrime =  phi - rod*kPhiSlicesPerSector*gridSizePhi ;
 772
 773
 774           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi
 775
 776           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
 777
 778             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 779                                               rlist, zedlist, philist, fLookUpBasic4ErOverEz  );
 780             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 781                                               rlist, zedlist, philist, fLookUpBasic4EphiOverEz);
 782             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 783                                               rlist, zedlist, philist, fLookUpBasic4DeltaEz   );
 784
 785           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
 786
 787             phiPrime     = TMath::TwoPi() - phiPrime ;
 788
 789             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 790                                               rlist, zedlist, philist, fLookUpBasic4ErOverEz  );
 791             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 792                                               rlist, zedlist, philist, fLookUpBasic4EphiOverEz);
 793             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 794                                               rlist, zedlist, philist, fLookUpBasic4DeltaEz   );
 795
 796             // Flip symmetry of phi integral for symmetric boundary conditions
 797             if ( symmetry[3] ==  1 ) ephiIntegral *= -1  ;
 798             // Flip symmetry of r integral if anti-symmetric boundary conditions
 799             if ( symmetry[3] == -1 ) erIntegral   *= -1  ;
 800
 801           }
 802
 803           if ( fgkZList[i] > 0 ) {
 804             erIntegralSum   += fRotatedClipVoltA[1]*erIntegral   ;
 805             ephiIntegralSum += fRotatedClipVoltA[1]*ephiIntegral ;
 806             ezIntegralSum   += fRotatedClipVoltA[1]*ezIntegral;
 807           } else if ( fgkZList[i] < 0 ) {
 808             erIntegralSum   += fRotatedClipVoltC[1]*erIntegral   ;
 809             ephiIntegralSum += fRotatedClipVoltC[1]*ephiIntegral ;
 810             ezIntegralSum   -= fRotatedClipVoltC[1]*ezIntegral;
 811           }
 812         }
 813
 814         // IFC Cooper rod shift  +++++++++++++++++++++++++++++
 815         for ( Int_t rod = 0 ; rod < 18 ; rod++ ) {
 816
 817           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
 818
 819           if ( fCopperRodShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
 820           if ( fCopperRodShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
 821
 822           // Rotate to a position where we have maps and prepare to sum
 823           phiPrime =  phi - rod*kPhiSlicesPerSector*gridSizePhi ;
 824
 825           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi
 826
 827           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
 828
 829             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 830                                               rlist, zedlist, philist, fLookUpBasic5ErOverEz  );
 831             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 832                                               rlist, zedlist, philist, fLookUpBasic5EphiOverEz);
 833             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 834                                               rlist, zedlist, philist, fLookUpBasic5DeltaEz   );
 835
 836           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
 837
 838             phiPrime     = TMath::TwoPi() - phiPrime ;
 839
 840             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 841                                               rlist, zedlist, philist, fLookUpBasic5ErOverEz  );
 842             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 843                                               rlist, zedlist, philist, fLookUpBasic5EphiOverEz);
 844             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 845                                               rlist, zedlist, philist, fLookUpBasic5DeltaEz   );
 846
 847             // Flip symmetry of phi integral for symmetric boundary conditions
 848             if ( symmetry[4] ==  1 ) ephiIntegral *= -1  ;
 849             // Flip symmetry of r integral if anti-symmetric boundary conditions
 850             if ( symmetry[4] == -1 ) erIntegral   *= -1  ;
 851
 852           }
 853
 854           if ( fgkZList[i] > 0 ) {
 855             erIntegralSum   += fCopperRodShiftA[rod]*erIntegral   ;
 856             ephiIntegralSum += fCopperRodShiftA[rod]*ephiIntegral ;
 857             ezIntegralSum   += fCopperRodShiftA[rod]*ezIntegral;
 858           } else if ( fgkZList[i] < 0 ) {
 859             erIntegralSum   += fCopperRodShiftC[rod]*erIntegral   ;
 860             ephiIntegralSum += fCopperRodShiftC[rod]*ephiIntegral ;
 861             ezIntegralSum   -= fCopperRodShiftC[rod]*ezIntegral;
 862           }
 863         }
 864
 865         // OFC Cooper rod shift  +++++++++++++++++++++++++++++
 866         for ( Int_t rod = 18 ; rod < 36 ; rod++ ) {
 867
 868           erIntegral = ephiIntegral = ezIntegral = 0.0 ;
 869
 870           if ( fCopperRodShiftA[rod] == 0 && fgkZList[i] > 0) continue ;
 871           if ( fCopperRodShiftC[rod] == 0 && fgkZList[i] < 0) continue ;
 872
 873           // Rotate to a position where we have maps and prepare to sum
 874           phiPrime =  phi - (rod-18)*kPhiSlicesPerSector*gridSizePhi ;
 875
 876           if ( phiPrime < 0 ) phiPrime += TMath::TwoPi() ;   // Stay in range from 0 to TwoPi
 877
 878           if ( (phiPrime >= 0) && (phiPrime <= kPhiSlices*gridSizePhi) ) {
 879
 880             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 881                                               rlist, zedlist, philist, fLookUpBasic6ErOverEz  );
 882             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 883                                               rlist, zedlist, philist, fLookUpBasic6EphiOverEz);
 884             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 885                                               rlist, zedlist, philist, fLookUpBasic6DeltaEz   );
 886
 887           }  else if ( (phiPrime <= TMath::TwoPi()) && (phiPrime >= (TMath::TwoPi()-kPhiSlices*gridSizePhi)) ){
 888
 889             phiPrime     = TMath::TwoPi() - phiPrime ;
 890
 891             erIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 892                                               rlist, zedlist, philist, fLookUpBasic6ErOverEz  );
 893             ephiIntegral = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 894                                               rlist, zedlist, philist, fLookUpBasic6EphiOverEz);
 895             ezIntegral   = Interpolate3DTable(order, r, z, phiPrime, kRows, kColumns, kPhiSlices,
 896                                               rlist, zedlist, philist, fLookUpBasic6DeltaEz   );
 897
 898             // Flip symmetry of phi integral for symmetric boundary conditions
 899             if ( symmetry[5] ==  1 ) ephiIntegral *= -1  ;
 900             // Flip symmetry of r integral if anti-symmetric boundary conditions
 901             if ( symmetry[5] == -1 ) erIntegral   *= -1  ;
 902
 903           }
 904
 905           if ( fgkZList[i] > 0 ) {
 906             erIntegralSum   += fCopperRodShiftA[rod]*erIntegral   ;
 907             ephiIntegralSum += fCopperRodShiftA[rod]*ephiIntegral ;
 908             ezIntegralSum   += fCopperRodShiftA[rod]*ezIntegral;
 909           } else if ( fgkZList[i] < 0 ) {
 910             erIntegralSum   += fCopperRodShiftC[rod]*erIntegral   ;
 911             ephiIntegralSum += fCopperRodShiftC[rod]*ephiIntegral ;
 912             ezIntegralSum   -= fCopperRodShiftC[rod]*ezIntegral;
 913           }
 914         }
 915
 916         // put the sum into the final lookup table
 917         erOverEz(j,i)   =  erIntegralSum;
 918         ephiOverEz(j,i) =  ephiIntegralSum;
 919         deltaEz(j,i)    =  ezIntegralSum;
 920
 921         //      if (j==1) printf("%lf %lf %lf: %lf %lf %lf\n",r, phi, z, erIntegralSum,ephiIntegralSum,ezIntegralSum);
 922
 923       } // endo r loop
 924     } // end of z loop
 925   } // end of phi loop
 926
 927
 928   // clear the temporary arrays lists
 929   for ( Int_t k = 0 ; k < kPhiSlices ; k++ )  {
 930     delete arrayofArrayV[k];
 931     delete arrayofCharge[k];
 932   }
 933
 934   AliInfo(" done");
 935   fInitLookUp = kTRUE;
 936
 937 }
 938
 939 void AliTPCFCVoltError3D::Print(const Option_t* option) const {
 940   //
 941   // Print function to check the settings of the Rod shifts and the rotated clips
 942   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 943   //
 944
 945   TString opt = option; opt.ToLower();
 946   printf("%s\n",GetTitle());
 947   printf(" - ROD shifts  (residual voltage settings): 40V correspond to 1mm of shift\n");
 948   Int_t count = 0;
 949   printf("  : A-side - (Rod & Strips) \n     ");
 950   for (Int_t i=0; i<36;i++) {
 951     if (fRodVoltShiftA[i]!=0) {
 952       printf("Rod%2d:%3.1fV ",i,fRodVoltShiftA[i]);
 953       count++;
 954     }
 955     if (count==0&&i==35)
 956       printf("-> all at 0 V\n");
 957     else if (i==35){
 958       printf("\n");
 959       count=0;
 960     }
 961   }
 962   printf("  : C-side - (Rod & Strips) \n     ");
 963   for (Int_t i=0; i<36;i++) {
 964     if (fRodVoltShiftC[i]!=0) {
 965       printf("Rod%2d:%3.1fV ",i,fRodVoltShiftC[i]);
 966       count++;
 967     }
 968     if (count==0&&i==35)
 969       printf("-> all at 0 V\n");
 970     else if (i==35){
 971       printf("\n");
 972       count=0;
 973     }
 974   }
 975
 976   printf(" - Rotated clips (residual voltage settings): 40V correspond to 1mm of 'offset'\n");
 977   if (fRotatedClipVoltA[0]!=0) {   printf("     A side (IFC): HV Rod: %3.1f V \n",fRotatedClipVoltA[0]); count++; }
 978   if (fRotatedClipVoltA[1]!=0) {   printf("     A side (OFC): HV Rod: %3.1f V \n",fRotatedClipVoltA[1]); count++; }
 979   if (fRotatedClipVoltC[0]!=0) {   printf("     C side (IFC): HV Rod: %3.1f V \n",fRotatedClipVoltC[0]); count++; }
 980   if (fRotatedClipVoltC[1]!=0) {   printf("     C side (OFC): HV Rod: %3.1f V \n",fRotatedClipVoltC[1]); count++; }
 981   if (count==0)
 982     printf("     -> no rotated clips \n");
 983
 984   count=0;
 985   printf(" - Copper ROD shifts (without strips):\n");
 986   printf("  : A-side - (Copper Rod shift) \n     ");
 987   for (Int_t i=0; i<36;i++) {
 988     if (fCopperRodShiftA[i]!=0) {
 989       printf("Rod%2d:%3.1fV ",i,fCopperRodShiftA[i]);
 990       count++;
 991     }
 992     if (count==0&&i==35)
 993       printf("-> all at 0 V\n");
 994     else if (i==35){
 995       printf("\n");
 996       count=0;
 997     }
 998   }
 999   printf("  : C-side - (Copper Rod shift) \n     ");
1000   for (Int_t i=0; i<36;i++) {
1001     if (fCopperRodShiftC[i]!=0) {
1002       printf("Rod%2d:%3.1fV ",i,fCopperRodShiftC[i]);
1003       count++;
1004     }
1005     if (count==0&&i==35)
1006       printf("-> all at 0 V\n");
1007     else if (i==35){
1008       printf("\n");
1009       count=0;
1010     }
1011   }
1012
1013   if (opt.Contains("a")) { // Print all details
1014     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
1015     printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
1016   }
1017
1018   if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitFCVoltError3D() ...");
1019
1020 }