TPC/TPCbase/AliTPCPointCorrection.cxx

   1 /// \class AliTPCPointCorrection
   2
   3 /**************************************************************************
   4  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   5  *                                                                        *
   6  * Author: The ALICE Off-line Project.                                    *
   7  * Contributors are mentioned in the code where appropriate.              *
   8  *                                                                        *
   9  * Permission to use, copy, modify and distribute this software and its   *
  10  * documentation strictly for non-commercial purposes is hereby granted   *
  11  * without fee, provided that the above copyright notice appears in all   *
  12  * copies and that both the copyright notice and this permission notice   *
  13  * appear in the supporting documentation. The authors make no claims     *
  14  * about the suitability of this software for any purpose. It is          *
  15  * provided "as is" without express or implied warranty.                  *
  16  **************************************************************************/
  17
  18 /*
  19
  20   Unlinearities fitting:
  21
  22   Model for Outer field cage:
  23   Unlinearities at the edge aproximated using two exponential decays.
  24
  25   dz = dz0(r,z) +dr(r,z)*tan(theta)
  26   dy =          +dr(r,z)*tan(phi)
  27
  28
  29
  30
  31 */
  32
  33 #include "AliTPCcalibDB.h"
  34 #include "TLinearFitter.h"
  35 #include "Riostream.h"
  36 #include "TList.h"
  37 #include "TMath.h"
  38 #include "TCanvas.h"
  39 #include "TFile.h"
  40 #include "TF1.h"
  41 #include "TVectorD.h"
  42 #include "AliLog.h"
  43 #include "AliTPCROC.h"
  44 #include "AliTPCClusterParam.h"
  45 #include "AliTPCPointCorrection.h"
  46
  47 AliTPCPointCorrection* AliTPCPointCorrection::fgInstance = 0;
  48
  49 /// \cond CLASSIMP
  50 ClassImp(AliTPCPointCorrection)
  51 /// \endcond
  52
  53 AliTPCPointCorrection::AliTPCPointCorrection():
  54   TNamed(),
  55   fParamsOutR(38),
  56   fParamsOutZ(38),
  57   fParamOutRVersion(0),
  58   fErrorsOutR(38),
  59   fErrorsOutZ(38),
  60   fParamOutZVersion(0),
  61   //
  62   fXIO(0),
  63   fXmiddle(0),
  64   fXquadrant(0),
  65   fArraySectorIntParam(36), // array of sector alignment parameters
  66   fArraySectorIntCovar(36), // array of sector alignment covariances
  67   //
  68   // Kalman filter for global alignment
  69   //
  70   fSectorParam(0),     // Kalman parameter
  71   fSectorCovar(0)     // Kalman covariance
  72
  73 {
  74   //
  75   // Default constructor
  76   //
  77   AliTPCROC * roc = AliTPCROC::Instance();
  78   fXquadrant = roc->GetPadRowRadii(36,53);
  79   fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
  80   fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
  81 }
  82
  83 AliTPCPointCorrection::AliTPCPointCorrection(const Text_t *name, const Text_t *title):
  84   TNamed(name,title),
  85   fParamsOutR(38),
  86   fParamsOutZ(38),
  87   fParamOutRVersion(0),
  88   fErrorsOutR(38),
  89   fErrorsOutZ(38),
  90   fParamOutZVersion(0),
  91   //
  92   //
  93   //
  94   fXIO(0),
  95   fXmiddle(0),
  96   fXquadrant(0),
  97   fArraySectorIntParam(36), // array of sector alignment parameters
  98   fArraySectorIntCovar(36), // array of sector alignment covariances
  99   //
 100   // Kalman filter for global alignment
 101   //
 102   fSectorParam(0),     // Kalman parameter   for A side
 103   fSectorCovar(0)     // Kalman covariance  for A side
 104
 105 {
 106   /// Non default constructor
 107
 108   AliTPCROC * roc = AliTPCROC::Instance();
 109   fXquadrant = roc->GetPadRowRadii(36,53);
 110   fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
 111   fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
 112
 113 }
 114
 115 AliTPCPointCorrection::~AliTPCPointCorrection(){
 116   ///
 117
 118 }
 119
 120
 121 AliTPCPointCorrection* AliTPCPointCorrection::Instance()
 122 {
 123   /// Singleton implementation
 124   /// Returns an instance of this class, it is created if neccessary
 125
 126   if (fgInstance == 0){
 127     fgInstance = new AliTPCPointCorrection();
 128   }
 129   return fgInstance;
 130 }
 131
 132
 133
 134 Double_t AliTPCPointCorrection::GetDrOut(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
 135   /// return radial correction
 136
 137   if (fParamOutRVersion==0) return CorrectionOutR0(isGlobal, type,cx,cy,cz,sector);
 138   return 0;
 139 }
 140
 141 Double_t      AliTPCPointCorrection::SGetDrOut(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
 142   /// return radial correction - static function
 143
 144   return fgInstance->GetDrOut(isGlobal, type,cx,cy,cz,sector);
 145 }
 146
 147
 148
 149
 150 Double_t AliTPCPointCorrection::GetDzOut(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
 151   ///
 152
 153   if (fParamOutZVersion==0) return CorrectionOutZ0(isGlobal, type,cx,cy,cz,sector);
 154   return 0;
 155 }
 156
 157
 158 Double_t      AliTPCPointCorrection::SGetDzOut(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector){
 159   ///
 160
 161   return fgInstance->GetDzOut(isGlobal, type,cx,cy,cz,sector);
 162 }
 163
 164
 165
 166
 167 Double_t AliTPCPointCorrection::CorrectionOutR0(Bool_t isGlobal, Bool_t type,  Double_t cx, Double_t cy, Double_t cz, Int_t sector){
 168   /// return dR correction - for correction version 0
 169   /// Parameters:
 170   /// isGlobal   - is the x in global frame
 171   /// type       - kTRUE - use sectors - kFALSE - only Side param
 172   /// cx, cy,cz  - cluster position
 173   /// sector     - sector number
 174
 175   if (isGlobal){
 176     // recalculate sector if in global frame
 177     Double_t alpha    = TMath::ATan2(cy,cx);
 178     if (alpha<0)  alpha+=TMath::Pi()*2;
 179     sector = Int_t(18*alpha/(TMath::Pi()*2));
 180   }
 181
 182   if (type==kFALSE) sector=36+(sector%36>=18);  // side level parameters
 183   // dout
 184   Double_t radius = TMath::Sqrt(cx*cx+cy*cy);
 185   AliTPCROC * roc = AliTPCROC::Instance();
 186   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
 187   Double_t dout = xout-radius;
 188   if (dout<0) return 0;
 189   //drift
 190   Double_t dr   = 0.5 - TMath::Abs(cz/250.);
 191   //
 192   //
 193   TVectorD * vec = GetParamOutR(sector);
 194   if (!vec) return 0;
 195   Double_t eout10 = TMath::Exp(-dout/10.);
 196   Double_t eout5 = TMath::Exp(-dout/5.);
 197   Double_t eoutp  = (eout10+eout5)*0.5;
 198   Double_t eoutm  = (eout10-eout5)*0.5;
 199   //
 200   Double_t result=0;
 201   result+=(*vec)[1]*eoutp;
 202   result+=(*vec)[2]*eoutm;
 203   result+=(*vec)[3]*eoutp*dr;
 204   result+=(*vec)[4]*eoutm*dr;
 205   result+=(*vec)[5]*eoutp*dr*dr;
 206   result+=(*vec)[6]*eoutm*dr*dr;
 207   return result;
 208 }
 209
 210 Double_t AliTPCPointCorrection::RPhiCOGCorrection(Int_t sector, Int_t padrow, Float_t pad, Float_t cy, Float_t y, Float_t z, Float_t ky,Float_t qMax, Float_t threshold){
 211   /// Calculates COG corection in RPHI direction
 212   /// cluster and track position  y is supposed to be corrected before for other effects
 213   /// (e.g ExB and alignemnt)
 214   /// Rphi distortion dependeds on the distance to the edge-pad, distance to the wire edge and
 215   /// relative distance to the center of the pad. Therefore the y position is trnsfromed to the
 216   /// pad coordiante frame (correction offset (ExB alignemnt) substracted).
 217   ///
 218   /// Input parameters:
 219   ///
 220   /// sector - sector number - 0-71  - cl.GetDetector()
 221   /// padrow - padrow number -0-63 -IROC 0-95 OROC - cl->GetRow()
 222   /// pad    - mean pad number  - cl->GetPad()
 223   /// cy     - cluster y        - cl->GetY()
 224   ///  y     - track -or cluster y
 225   ///  qMax  - cluster max charge - cl-.GetMax()
 226   ///  threshold - clusterer threshold
 227
 228   AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
 229   if (!clparam) {
 230     AliFatal("TPC OCDB not initialized");
 231     return 0;
 232   }
 233   Int_t padtype=0;
 234   if (sector>=36) padtype = (padrow>64)?2:1;
 235   Double_t padwidth=(padtype==0)? 0.4:0.6;
 236
 237   Double_t sigma = clparam->GetRMS0(0,padtype,250-TMath::Abs(z),ky)/padwidth;
 238   //
 239   Int_t   npads   =  AliTPCROC::Instance()->GetNPads(sector,padrow);
 240   Float_t yshift  =  TMath::Abs(cy)-TMath::Abs((-npads*0.5+pad)*padwidth);   // the clusters can be corrected before
 241                                             // shift to undo correction
 242
 243   y -= yshift*((y>0)?1.:-1.);                             // y in the sector coordinate
 244   Double_t y0     = npads*0.5*padwidth;
 245   Double_t dy     = (TMath::Abs(y0)-TMath::Abs(y))/padwidth-0.5;  // distance to  the center of pad0
 246   //
 247   Double_t padcenter = TMath::Nint(dy);
 248   Double_t sumw=0;
 249   Double_t sumwi=0;
 250   for (Double_t ip=padcenter-2.5; ip<=padcenter+2.5;ip++){
 251     Double_t weightGaus = qMax*TMath::Exp(-(dy-ip)*(dy-ip)/(2*sigma*sigma));
 252     Double_t ypad       = (ip-npads*0.5)*padwidth;
 253     Double_t weightGain = GetEdgeQ0(sector,padrow,ypad);
 254     //
 255     Double_t weight = TMath::Nint(weightGaus*weightGain);
 256     if (ip<0 &&weight> threshold) weight = threshold;  // this is done in cl finder
 257     if (weight<0) continue;
 258     sumw+=weight;
 259     sumwi+=weight*(ip);
 260   }
 261   Double_t result =0;
 262   if (sumw>0) result = sumwi/sumw;
 263   result = (result-dy)*padwidth;
 264   return result;
 265 }
 266
 267
 268
 269
 270 Double_t AliTPCPointCorrection::CorrectionOutZ0(Bool_t isGlobal, Bool_t type,  Double_t cx, Double_t cy, Double_t cz, Int_t sector){
 271   /// return dR correction - for correction version 0
 272   /// Parameters:
 273   /// isGlobal   - is the x in global frame
 274   /// type       - kTRUE - use sectors - kFALSE - only Side param
 275   /// cx, cy,cz  - cluster position
 276   /// sector     - sector number
 277
 278   if (isGlobal){
 279     // recalculate sector if in global frame
 280     Double_t alpha    = TMath::ATan2(cy,cx);
 281     if (alpha<0)  alpha+=TMath::Pi()*2;
 282     sector = Int_t(18*alpha/(TMath::Pi()*2));
 283   }
 284
 285   if (type==kFALSE) sector=36+(sector%36>=18);  // side level parameters
 286   // dout
 287   Double_t radius = TMath::Sqrt(cx*cx+cy*cy);
 288   AliTPCROC * roc = AliTPCROC::Instance();
 289   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
 290   Double_t dout = xout-radius;
 291   if (dout<0) return 0;
 292   //drift
 293   Double_t dr   = 0.5 - TMath::Abs(cz/250.);
 294   //
 295   //
 296   TVectorD * vec = GetParamOutR(sector);
 297   if (!vec) return 0;
 298   Double_t eout10 = TMath::Exp(-dout/10.);
 299   Double_t eout5 = TMath::Exp(-dout/5.);
 300   Double_t eoutp  = (eout10+eout5)*0.5;
 301   Double_t eoutm  = (eout10-eout5)*0.5;
 302   //
 303   Double_t result=0;
 304   result+=(*vec)[1]*eoutp;
 305   result+=(*vec)[2]*eoutm;
 306   result+=(*vec)[3]*eoutp*dr;
 307   result+=(*vec)[4]*eoutm*dr;
 308   result+=(*vec)[5]*eoutp*dr*dr;
 309   result+=(*vec)[6]*eoutm*dr*dr;
 310   return result;
 311
 312 }
 313
 314 Double_t  AliTPCPointCorrection::GetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
 315   ///
 316
 317   /* TF1 fexp("fexp","1-exp(-[0]*(x-[1]))",0,20)
 318           | param [0] | param [1]
 319      IROC | 4.71413   | 1.39558
 320      OROC1| 2.11437   | 1.52643
 321      OROC2| 2.15082   | 1.53537
 322   */
 323
 324   Double_t params[2]={100,0};
 325   if (sector<36){
 326     params[0]=4.71413; params[1]=1.39558;
 327   }
 328   if (sector>=36){
 329     if (padrow<64) {  params[0]=2.114; params[1]=1.526;}
 330     if (padrow>=64){  params[0]=2.15; params[1]=1.535;}
 331   }
 332   Double_t result= 1;
 333   Double_t xrow  = AliTPCROC::Instance()->GetPadRowRadii(sector,padrow);
 334   Double_t ymax  = TMath::Tan(TMath::Pi()/18.)*xrow;
 335   Double_t dedge = ymax-TMath::Abs(y);
 336   if (dedge-params[1]<0)             return 0;
 337   if (dedge>10.*params[1]) return 1;
 338   result = 1.-TMath::Exp(-params[0]*(dedge-params[1]));
 339   return result;
 340 }
 341
 342 Double_t AliTPCPointCorrection::SRPhiCOGCorrection(Int_t sector, Int_t padrow, Float_t pad, Float_t cy, Float_t y, Float_t z, Float_t ky,Float_t qMax, Float_t threshold){
 343   return fgInstance->RPhiCOGCorrection(sector, padrow, pad, cy, y,z, ky, qMax, threshold);
 344 }
 345
 346 Double_t AliTPCPointCorrection::SGetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
 347   ///
 348
 349   return fgInstance->GetEdgeQ0(sector, padrow, y);
 350 }
 351
 352 void     AliTPCPointCorrection::AddCorrectionSector(TObjArray & sideAPar, TObjArray &sideCPar, TObjArray & sideACov, TObjArray &sideCCov, Bool_t reset){
 353   ///
 354
 355   for (Int_t isec=0;isec<36;isec++){
 356     TMatrixD * mat1 = (TMatrixD*)(sideAPar.At(isec));
 357     TMatrixD * mat1Covar = (TMatrixD*)(sideACov.At(isec));
 358     if (!mat1) continue;
 359     TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
 360     TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
 361     if (!mat0) {
 362       fArraySectorIntParam.AddAt(mat1->Clone(),isec);
 363       fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec);
 364       continue;
 365     }
 366     (*mat0Covar)=(*mat1Covar);
 367     if (reset)   (*mat0)=(*mat1);
 368     if (!reset) (*mat0)+=(*mat1);
 369   }
 370
 371   for (Int_t isec=0;isec<36;isec++){
 372     TMatrixD * mat1 = (TMatrixD*)(sideCPar.At(isec));
 373     TMatrixD * mat1Covar = (TMatrixD*)(sideCCov.At(isec));
 374     if (!mat1) continue;
 375     TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
 376     TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
 377     if (!mat0) {
 378       fArraySectorIntParam.AddAt(mat1->Clone(),isec);
 379       fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec);
 380       continue;
 381     }
 382     (*mat0Covar)=(*mat1Covar);
 383     if (reset)   (*mat0)=(*mat1);
 384     if (!reset) (*mat0)+=(*mat1);
 385   }
 386 }
 387
 388
 389 Double_t AliTPCPointCorrection::GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/, Int_t quadrant){
 390   /// Get position correction for given sector
 391
 392   TMatrixD * param = (TMatrixD*)fArraySectorIntParam.At(sector%36);
 393   if (!param) return 0;
 394   if (quadrant<0){   //recalc quadrant if not specified
 395     if (lx<fXIO) quadrant=0;
 396     if(lx>fXIO){
 397       if (lx<fXquadrant) {
 398         if (ly<0) quadrant=1;
 399         if (ly>0) quadrant=2;
 400       }
 401       if (lx>=fXquadrant) {
 402         if (ly<0) quadrant=3;
 403         if (ly>0) quadrant=4;
 404       }
 405     }
 406   }
 407   Double_t a10 = (*param)(quadrant*6+0,0);
 408   Double_t a20 = (*param)(quadrant*6+1,0);
 409   Double_t a21 = (*param)(quadrant*6+2,0);
 410   Double_t dx  = (*param)(quadrant*6+3,0);
 411   Double_t dy  = (*param)(quadrant*6+4,0);
 412   Double_t dz  = (*param)(quadrant*6+5,0);
 413   Double_t deltaX = lx-fXmiddle;
 414   if (coord==0) return dx;
 415   if (coord==1) return dy+deltaX*a10;
 416   if (coord==2) return dz+deltaX*a20+ly*a21;
 417   if (coord==3) return a10;
 418   if (coord==4) return a20;
 419   if (coord==5) return a21;
 420   //
 421   return 0;
 422 }
 423
 424 Double_t AliTPCPointCorrection::SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz, Int_t quadrant){
 425   ///
 426
 427   if (!Instance()) return 0;
 428   return Instance()->GetCorrectionSector(coord,sector,lx,ly,lz, quadrant);
 429 }
 430
 431
 432
 433 Double_t AliTPCPointCorrection::GetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/){
 434   /// Get position correction for given sector
 435
 436   if (!fSectorParam) return 0;
 437
 438   Double_t a10 = (*fSectorParam)(sector*6+0,0);
 439   Double_t a20 = (*fSectorParam)(sector*6+1,0);
 440   Double_t a21 = (*fSectorParam)(sector*6+2,0);
 441   Double_t dx  = (*fSectorParam)(sector*6+3,0);
 442   Double_t dy  = (*fSectorParam)(sector*6+4,0);
 443   Double_t dz  = (*fSectorParam)(sector*6+5,0);
 444   Double_t deltaX = lx-fXIO;
 445   //
 446   if (coord==0) return dx;
 447   if (coord==1) return dy+deltaX*a10;
 448   if (coord==2) return dz+deltaX*a20+ly*a21;
 449   if (coord==3) return a10;
 450   if (coord==4) return a20;
 451   if (coord==5) return a21;
 452   return 0;
 453 }
 454
 455 Double_t AliTPCPointCorrection::SGetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){
 456   ///
 457
 458   if (!Instance()) return 0;
 459   return Instance()->GetCorrection(coord,sector,lx,ly,lz);
 460 }
 461
 462
 463
 464
 465
 466
 467