TPC/AliTPCROCVoltError3D.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 // AliTPCROCVoltError3D class                                               //
  19 // The class calculates the space point distortions due to residual voltage //
  20 // errors on Read Out Chambers 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 //       misalignment and the vd dependency on the residual drift field     //
  26 //                                                                          //
  27 // date: 08/08/2010                                                         //
  28 // Authors: Jim Thomas, Stefan Rossegger                                    //
  29 //                                                                          //
  30 // Example usage :                                                          //
  31 //  AliTPCROCVoltError3D ROCerror;                                            //
  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 #include "TFile.h"
  41
  42 #include "TMath.h"
  43 #include "AliTPCROC.h"
  44 #include "AliTPCROCVoltError3D.h"
  45
  46 ClassImp(AliTPCROCVoltError3D)
  47
  48 AliTPCROCVoltError3D::AliTPCROCVoltError3D()
  49   : AliTPCCorrection("ROCVoltErrors","ROC z alignment Errors"),
  50     fC0(0.),fC1(0.),
  51     fROCdisplacement(kTRUE),
  52     fInitLookUp(kFALSE),
  53     fROCDataFileName("$(ALICE_ROOT)/TPC/Calib/maps/TPCROCdzSurvey.root"),  // standard file name of ROC survey
  54     fdzDataLinFit(0)
  55 {
  56   //
  57   // default constructor
  58   //
  59
  60   // Array which will contain the solution according to the setted boundary conditions
  61   // main input: z alignment of the Read Out chambers
  62   // see InitROCVoltError3D() function
  63   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
  64     fLookUpErOverEz[k]   =  new TMatrixD(kNR,kNZ);
  65     fLookUpEphiOverEz[k] =  new TMatrixD(kNR,kNZ);
  66     fLookUpDeltaEz[k]    =  new TMatrixD(kNR,kNZ);
  67   }
  68
  69   SetROCDataFileName(fROCDataFileName); // initialization of fdzDataLinFit is included
  70
  71 }
  72
  73 AliTPCROCVoltError3D::~AliTPCROCVoltError3D() {
  74   //
  75   // destructor
  76   //
  77
  78   for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
  79     delete fLookUpErOverEz[k];
  80     delete fLookUpEphiOverEz[k];
  81     delete fLookUpDeltaEz[k];
  82   }
  83
  84   delete fdzDataLinFit;
  85 }
  86
  87 void AliTPCROCVoltError3D::Init() {
  88   //
  89   // Initialization funtion
  90   //
  91
  92   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  93   if (!magF) AliError("Magneticd field - not initialized");
  94   Double_t bzField = magF->SolenoidField()/10.; //field in T
  95   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
  96   if (!param) AliError("Parameters - not initialized");
  97   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
  98   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
  99   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 100   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 101   SetOmegaTauT1T2(wt,fT1,fT2);
 102
 103   InitROCVoltError3D();
 104 }
 105
 106 void AliTPCROCVoltError3D::Update(const TTimeStamp &/*timeStamp*/) {
 107   //
 108   // Update function
 109   //
 110   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 111   if (!magF) AliError("Magneticd field - not initialized");
 112   Double_t bzField = magF->SolenoidField()/10.; //field in T
 113   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 114   if (!param) AliError("Parameters - not initialized");
 115   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]  // From dataBase: to be updated: per second (ideally)
 116   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 117   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 118   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 119   SetOmegaTauT1T2(wt,fT1,fT2);
 120
 121 }
 122
 123 void  AliTPCROCVoltError3D::SetROCDataFileName(char *const fname) {
 124   //
 125   // Set / load the ROC data (linear fit of ROC misalignments)
 126   //
 127
 128   fROCDataFileName = fname;
 129
 130   TFile f(fROCDataFileName,"READ");
 131   TMatrixD *m = (TMatrixD*) f.Get("dzSurveyLinFitData");
 132   TMatrixD &mf = *m;
 133
 134   // prepare some space
 135
 136   if (fdzDataLinFit) delete fdzDataLinFit;
 137   fdzDataLinFit = new TMatrixD(72,3);
 138   TMatrixD &dataIntern = *fdzDataLinFit;
 139
 140   for (Int_t iroc=0;iroc<72;iroc++) {
 141     dataIntern(iroc,0) = mf(iroc,0);  // z0 offset
 142     dataIntern(iroc,1) = mf(iroc,1);  // slope in x
 143     dataIntern(iroc,2) = mf(iroc,2);  // slope in y
 144   }
 145
 146   f.Close();
 147
 148   fInitLookUp = kFALSE;
 149
 150 }
 151
 152 void AliTPCROCVoltError3D::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
 153   //
 154   // Calculates the correction due e.g. residual voltage errors on the TPC boundaries
 155   //
 156
 157   if (!fInitLookUp) {
 158     AliInfo("Lookup table was not initialized! Perform the inizialisation now ...");
 159     InitROCVoltError3D();
 160     return;
 161   }
 162
 163   Int_t   order     = 1 ;               // FIXME: hardcoded? Linear interpolation = 1, Quadratic = 2
 164
 165   Double_t intEr, intEphi, intDeltaEz;
 166   Double_t r, phi, z ;
 167   Int_t    sign;
 168
 169   r      =  TMath::Sqrt( x[0]*x[0] + x[1]*x[1] ) ;
 170   phi    =  TMath::ATan2(x[1],x[0]) ;
 171   if ( phi < 0 ) phi += TMath::TwoPi() ;                   // Table uses phi from 0 to 2*Pi
 172   z      =  x[2] ;                                         // Create temporary copy of x[2]
 173
 174   if ( (roc%36) < 18 ) {
 175     sign =  1;       // (TPC A side)
 176   } else {
 177     sign = -1;       // (TPC C side)
 178   }
 179
 180   if ( sign==1  && z <  fgkZOffSet ) z =  fgkZOffSet;    // Protect against discontinuity at CE
 181   if ( sign==-1 && z > -fgkZOffSet ) z = -fgkZOffSet;    // Protect against discontinuity at CE
 182
 183
 184   if ( (sign==1 && z<0) || (sign==-1 && z>0) ) // just a consistency check
 185     AliError("ROC number does not correspond to z coordinate! Calculation of distortions is most likely wrong!");
 186
 187   // Get the Er and Ephi field integrals plus the integral over DeltaEz
 188   intEr      = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 189                                   fgkRList, fgkZList, fgkPhiList, fLookUpErOverEz  );
 190   intEphi    = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 191                                   fgkRList, fgkZList, fgkPhiList, fLookUpEphiOverEz);
 192   intDeltaEz = Interpolate3DTable(order, r, z, phi, kNR, kNZ, kNPhi,
 193                                   fgkRList, fgkZList, fgkPhiList, fLookUpDeltaEz   );
 194
 195   //  printf("%lf %lf %lf\n",intEr,intEphi,intDeltaEz);
 196
 197   // Calculate distorted position
 198   if ( r > 0.0 ) {
 199     phi =  phi + ( fC0*intEphi - fC1*intEr ) / r;
 200     r   =  r   + ( fC0*intEr   + fC1*intEphi );
 201   }
 202
 203   // Calculate correction in cartesian coordinates
 204   dx[0] = r * TMath::Cos(phi) - x[0];
 205   dx[1] = r * TMath::Sin(phi) - x[1];
 206   dx[2] = intDeltaEz;  // z distortion - (internally scaled with driftvelocity dependency
 207                        // on the Ez field plus the actual ROC misalignment (if set TRUE)
 208
 209 }
 210
 211 void AliTPCROCVoltError3D::InitROCVoltError3D() {
 212   //
 213   // Initialization of the Lookup table which contains the solutions of the
 214   // Dirichlet boundary problem
 215   // Calculation of the single 3D-Poisson solver is done just if needed
 216   // (see basic lookup tables in header file)
 217   //
 218
 219   const Int_t   order       =    1  ;  // Linear interpolation = 1, Quadratic = 2
 220   const Float_t gridSizeR   =  (fgkOFCRadius-fgkIFCRadius) / (kRows-1) ;
 221   const Float_t gridSizeZ   =  fgkTPCZ0 / (kColumns-1) ;
 222   const Float_t gridSizePhi =  TMath::TwoPi() / ( 18.0 * kPhiSlicesPerSector);
 223
 224   // temporary arrays to create the boundary conditions
 225   TMatrixD *arrayofArrayV[kPhiSlices], *arrayofCharge[kPhiSlices] ;
 226   TMatrixD *arrayofEroverEz[kPhiSlices], *arrayofEphioverEz[kPhiSlices], *arrayofDeltaEz[kPhiSlices] ;
 227
 228   for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
 229     arrayofArrayV[k]     =   new TMatrixD(kRows,kColumns) ;
 230     arrayofCharge[k]     =   new TMatrixD(kRows,kColumns) ;
 231     arrayofEroverEz[k]   =   new TMatrixD(kRows,kColumns) ;
 232     arrayofEphioverEz[k] =   new TMatrixD(kRows,kColumns) ;
 233     arrayofDeltaEz[k]    =   new TMatrixD(kRows,kColumns) ;
 234   }
 235
 236   // list of point as used in the poisson relation and the interpolation (during sum up)
 237   Double_t  rlist[kRows], zedlist[kColumns] , philist[kPhiSlices];
 238   for ( Int_t k = 0 ; k < kPhiSlices ; k++ ) {
 239     philist[k] =  gridSizePhi * k;
 240     for ( Int_t i = 0 ; i < kRows ; i++ )    {
 241       rlist[i] = fgkIFCRadius + i*gridSizeR ;
 242       for ( Int_t j = 0 ; j < kColumns ; j++ ) { // Fill Vmatrix with Boundary Conditions
 243         zedlist[j]  = j * gridSizeZ ;
 244       }
 245     }
 246   }
 247
 248   // ==========================================================================
 249   // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
 250   // Allow for different size grid spacing in R and Z directions
 251
 252   const Int_t   symmetry = 0;
 253
 254   // Set bondaries and solve Poisson's equation --------------------------
 255
 256   if ( !fInitLookUp ) {
 257
 258     AliInfo(Form("Solving the poisson equation (~ %d sec)",2*10*(int)(kPhiSlices/10)));
 259
 260     for ( Int_t side = 0 ; side < 2 ; side++ ) {  // Solve Poisson3D twice; once for +Z and once for -Z
 261
 262       for ( Int_t k = 0 ; k < kPhiSlices ; k++ )  {
 263         TMatrixD &arrayV    =  *arrayofArrayV[k] ;
 264         TMatrixD &charge    =  *arrayofCharge[k] ;
 265
 266         //Fill arrays with initial conditions.  V on the boundary and Charge in the volume.
 267         for ( Int_t i = 0 ; i < kRows ; i++ ) {
 268           for ( Int_t j = 0 ; j < kColumns ; j++ ) {  // Fill Vmatrix with Boundary Conditions
 269             arrayV(i,j) = 0.0 ;
 270             charge(i,j) = 0.0 ;
 271
 272             Float_t radius0 = rlist[i] ;
 273             Float_t phi0    = gridSizePhi * k ;
 274
 275             // To avoid problems at sector boundaries, use an average of +- 1 degree from actual phi location
 276             if ( j == (kColumns-1) )
 277               arrayV(i,j) = 0.5*  ( GetROCVoltOffset( side, radius0, phi0+0.02 ) + GetROCVoltOffset( side, radius0, phi0-0.02 ) ) ;
 278
 279           }
 280         }
 281
 282         for ( Int_t i = 1 ; i < kRows-1 ; i++ ) {
 283           for ( Int_t j = 1 ; j < kColumns-1 ; j++ ) {
 284             charge(i,j)  =  0.0 ;
 285           }
 286         }
 287       }
 288
 289       // Solve Poisson's equation in 3D cylindrical coordinates by relaxation technique
 290       // Allow for different size grid spacing in R and Z directions
 291
 292       PoissonRelaxation3D( arrayofArrayV, arrayofCharge,
 293                            arrayofEroverEz, arrayofEphioverEz, arrayofDeltaEz,
 294                            kRows, kColumns, kPhiSlices, gridSizePhi, kIterations,
 295                            symmetry, fROCdisplacement) ;
 296
 297
 298       //Interpolate results onto a custom grid which is used just for these calculations.
 299       Double_t  r, phi, z ;
 300       for ( Int_t k = 0 ; k < kNPhi ; k++ ) {
 301         phi = fgkPhiList[k] ;
 302
 303         TMatrixD &erOverEz   =  *fLookUpErOverEz[k]  ;
 304         TMatrixD &ephiOverEz =  *fLookUpEphiOverEz[k];
 305         TMatrixD &deltaEz    =  *fLookUpDeltaEz[k]   ;
 306
 307         for ( Int_t j = 0 ; j < kNZ ; j++ ) {
 308
 309           z = TMath::Abs(fgkZList[j]) ;  // Symmetric solution in Z that depends only on ABS(Z)
 310
 311           if ( side == 0 &&  fgkZList[j] < 0 ) continue; // Skip rest of this loop if on the wrong side
 312           if ( side == 1 &&  fgkZList[j] > 0 ) continue; // Skip rest of this loop if on the wrong side
 313
 314           for ( Int_t i = 0 ; i < kNR ; i++ ) {
 315             r = fgkRList[i] ;
 316
 317             // Interpolate basicLookup tables; once for each rod, then sum the results
 318             erOverEz(i,j)   = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
 319                                                  rlist, zedlist, philist, arrayofEroverEz  );
 320             ephiOverEz(i,j) = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
 321                                                  rlist, zedlist, philist, arrayofEphioverEz);
 322             deltaEz(i,j)    = Interpolate3DTable(order, r, z, phi, kRows, kColumns, kPhiSlices,
 323                                                  rlist, zedlist, philist, arrayofDeltaEz  );
 324
 325             if (side == 1)  deltaEz(i,j) = -  deltaEz(i,j); // negative coordinate system on C side
 326
 327           } // end r loop
 328         }// end z loop
 329       }// end phi loop
 330
 331       if ( side == 0 ) AliInfo(" A side done");
 332       if ( side == 1 ) AliInfo(" C side done");
 333     } // end side loop
 334   }
 335
 336   // clear the temporary arrays lists
 337   for ( Int_t k = 0 ; k < kPhiSlices ; k++ )  {
 338     delete arrayofArrayV[k];
 339     delete arrayofCharge[k];
 340     delete arrayofEroverEz[k];
 341     delete arrayofEphioverEz[k];
 342     delete arrayofDeltaEz[k];
 343   }
 344
 345
 346   fInitLookUp = kTRUE;
 347
 348 }
 349
 350
 351 Float_t AliTPCROCVoltError3D::GetROCVoltOffset(Int_t side, Float_t r0, Float_t phi0) {
 352   //
 353   // Returns the dz alignment data (in voltage equivalents) at
 354   // the given position
 355   //
 356
 357   Float_t xp = r0*TMath::Cos(phi0);
 358   Float_t yp = r0*TMath::Sin(phi0);
 359
 360   // phi0 should be between 0 and 2pi
 361   if (phi0<0) phi0+=TMath::TwoPi();
 362   Int_t roc = (Int_t)TMath::Floor((TMath::RadToDeg()*phi0)/20);
 363   if (side==1) roc+=18; // C side
 364   if (r0>132) roc+=36;  // OROC
 365
 366   // linear-plane data:  z = z0 + kx*x + ky*y
 367   TMatrixD &fitData = *fdzDataLinFit;
 368   Float_t dz = fitData(roc,0)+fitData(roc,1)*xp + fitData(roc,2)*yp; // value in cm
 369
 370   // aproximated Voltage-offset-aquivalent to the z misalignment
 371   // (linearly scaled with the z position)
 372   Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
 373   Float_t voltOff = dz*ezField;            // values in "Volt equivalents"
 374
 375   return voltOff;
 376 }
 377
 378 TH2F * AliTPCROCVoltError3D::CreateHistoOfZSurvey(Int_t side, Int_t nx, Int_t ny) {
 379   //
 380   // return a simple histogramm containing the input to the poisson solver
 381   // (z positions of the Read-out chambers, linearly interpolated)
 382
 383   char hname[100];
 384   if (side==0) sprintf(hname,"survey_dz_Aside");
 385   if (side==1) sprintf(hname,"survey_dz_Cside");
 386
 387   TH2F *h = new TH2F(hname,hname,nx,-250.,250.,ny,-250.,250.);
 388
 389   for (Int_t iy=1;iy<=ny;++iy) {
 390     Double_t yp = h->GetYaxis()->GetBinCenter(iy);
 391     for (Int_t ix=1;ix<=nx;++ix) {
 392       Double_t xp = h->GetXaxis()->GetBinCenter(ix);
 393
 394       Float_t r0 = TMath::Sqrt(xp*xp+yp*yp);
 395       Float_t phi0 = TMath::ATan2(yp,xp);
 396
 397       Float_t dz = GetROCVoltOffset(side,r0,phi0); // in [volt]
 398
 399       Double_t ezField = (fgkCathodeV-fgkGG)/fgkTPCZ0; // = ALICE Electric Field (V/cm) Magnitude ~ -400 V/cm;
 400       dz = dz/ezField;    // in [cm]
 401
 402       if (85.<=r0 && r0<=245.) {
 403         h->SetBinContent(ix,iy,dz);
 404       } else {
 405         h->SetBinContent(ix,iy,0.);
 406       }
 407     }
 408   }
 409
 410   h->GetXaxis()->SetTitle("x [cm]");
 411   h->GetYaxis()->SetTitle("y [cm]");
 412   h->GetZaxis()->SetTitle("dz [cm]");
 413   h->SetStats(0);
 414   //  h->DrawCopy("colz");
 415
 416   return h;
 417 }
 418
 419 void AliTPCROCVoltError3D::Print(const Option_t* option) const {
 420   //
 421   // Print function to check the settings of the Rod shifts and the rotated clips
 422   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 423   //
 424
 425   TString opt = option; opt.ToLower();
 426   printf("%s\n",GetTitle());
 427   printf(" - Voltage settings on the TPC Read-Out chambers - linearly interpolated\n");
 428   printf("   info: Check the following data-file for more details: %s \n",fROCDataFileName);
 429
 430   if (opt.Contains("a")) { // Print all details
 431     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
 432     printf(" - C1: %1.4f, C0: %1.4f \n",fC1,fC0);
 433   }
 434
 435   if (!fInitLookUp) AliError("Lookup table was not initialized! You should do InitROCVoltError3D() ...");
 436
 437 }