TPC/AliTPCFCVoltError3D.cxx

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