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