TPC/Rec/AliTPCseed.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
  19 //-----------------------------------------------------------------
  20 //
  21 //           Implementation of the TPC seed class
  22 //        This class is used by the AliTPCtracker class
  23 //      Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
  24 //-----------------------------------------------------------------
  25 #include <TVectorF.h>
  26 #include "TClonesArray.h"
  27 #include "TGraphErrors.h"
  28 #include "AliTPCseed.h"
  29 #include "AliTPCReconstructor.h"
  30 #include "AliTPCClusterParam.h"
  31 #include "AliTPCCalPad.h"
  32 #include "AliTPCCalROC.h"
  33 #include "AliTPCcalibDB.h"
  34 #include "AliTPCParam.h"
  35 #include "AliMathBase.h"
  36 #include "AliTPCTransform.h"
  37 #include "AliSplineFit.h"
  38 #include "AliCDBManager.h"
  39 #include "AliTPCcalibDButil.h"
  40
  41
  42 ClassImp(AliTPCseed)
  43
  44
  45
  46 AliTPCseed::AliTPCseed():
  47   AliTPCtrack(),
  48   fEsd(0x0),
  49   fClusterOwner(kFALSE),
  50   fRow(0),
  51   fSector(-1),
  52   fRelativeSector(-1),
  53   fCurrentSigmaY2(1e10),
  54   fCurrentSigmaZ2(1e10),
  55   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
  56   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
  57   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
  58   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
  59   //
  60   fErrorY2(1e10),
  61   fErrorZ2(1e10),
  62   fCurrentCluster(0x0),
  63   fCurrentClusterIndex1(-1),
  64   fInDead(kFALSE),
  65   fIsSeeding(kFALSE),
  66   fNoCluster(0),
  67   fSort(0),
  68   fBSigned(kFALSE),
  69   fSeedType(0),
  70   fSeed1(-1),
  71   fSeed2(-1),
  72   fMAngular(0),
  73   fCircular(0),
  74   fPoolID(-1)
  75 {
  76   //
  77   for (Int_t i=0;i<160;i++) SetClusterIndex2(i,-3);
  78   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
  79   for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
  80   for (Int_t i=0;i<AliPID::kSPECIES;i++)   fTPCr[i]=0.2;
  81   for (Int_t i=0;i<4;i++) {
  82     fDEDX[i] = 0.;
  83     fSDEDX[i] = 1e10;
  84     fNCDEDX[i] = 0;
  85     fNCDEDXInclThres[i] = 0;
  86   }
  87   fDEDX[4] = 0;
  88   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
  89 }
  90
  91 AliTPCseed::AliTPCseed(const AliTPCseed &s, Bool_t clusterOwner):
  92   AliTPCtrack(s),
  93   fEsd(0x0),
  94   fClusterOwner(clusterOwner),
  95   fRow(0),
  96   fSector(-1),
  97   fRelativeSector(-1),
  98   fCurrentSigmaY2(-1),
  99   fCurrentSigmaZ2(-1),
 100   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
 101   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
 102   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
 103   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
 104   fErrorY2(1e10),
 105   fErrorZ2(1e10),
 106   fCurrentCluster(0x0),
 107   fCurrentClusterIndex1(-1),
 108   fInDead(kFALSE),
 109   fIsSeeding(kFALSE),
 110   fNoCluster(0),
 111   fSort(0),
 112   fBSigned(kFALSE),
 113   fSeedType(0),
 114   fSeed1(-1),
 115   fSeed2(-1),
 116   fMAngular(0),
 117   fCircular(0),
 118   fPoolID(-1)
 119 {
 120   //---------------------
 121   // dummy copy constructor
 122   //-------------------------
 123   for (Int_t i=0;i<160;i++) {
 124     fClusterPointer[i]=0;
 125     if (fClusterOwner){
 126       if (s.fClusterPointer[i])
 127         fClusterPointer[i] = new AliTPCclusterMI(*(s.fClusterPointer[i]));
 128     }else{
 129       fClusterPointer[i] = s.fClusterPointer[i];
 130     }
 131     fTrackPoints[i] = s.fTrackPoints[i];
 132   }
 133   for (Int_t i=0;i<160;i++) fIndex[i] = s.fIndex[i];
 134   for (Int_t i=0;i<AliPID::kSPECIES;i++)   fTPCr[i]=s.fTPCr[i];
 135   for (Int_t i=0;i<4;i++) {
 136     fDEDX[i] = s.fDEDX[i];
 137     fSDEDX[i] = s.fSDEDX[i];
 138     fNCDEDX[i] = s.fNCDEDX[i];
 139     fNCDEDXInclThres[i] = s.fNCDEDXInclThres[i];
 140   }
 141   fDEDX[4] = s.fDEDX[4];
 142   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = s.fOverlapLabels[i];
 143
 144 }
 145
 146
 147 AliTPCseed::AliTPCseed(const AliTPCtrack &t):
 148   AliTPCtrack(t),
 149   fEsd(0x0),
 150   fClusterOwner(kFALSE),
 151   fRow(0),
 152   fSector(-1),
 153   fRelativeSector(-1),
 154   fCurrentSigmaY2(-1),
 155   fCurrentSigmaZ2(-1),
 156   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
 157   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
 158   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
 159   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
 160   fErrorY2(1e10),
 161   fErrorZ2(1e10),
 162   fCurrentCluster(0x0),
 163   fCurrentClusterIndex1(-1),
 164   fInDead(kFALSE),
 165   fIsSeeding(kFALSE),
 166   fNoCluster(0),
 167   fSort(0),
 168   fBSigned(kFALSE),
 169   fSeedType(0),
 170   fSeed1(-1),
 171   fSeed2(-1),
 172   fMAngular(0),
 173   fCircular(0),
 174   fPoolID(-1)
 175 {
 176   //
 177   // Constructor from AliTPCtrack
 178   //
 179   fFirstPoint =0;
 180   for (Int_t i=0;i<5;i++)   fTPCr[i]=0.2;
 181   for (Int_t i=0;i<160;i++) {
 182     fClusterPointer[i] = 0;
 183     Int_t index = t.GetClusterIndex(i);
 184     if (index>=-1){
 185       SetClusterIndex2(i,index);
 186     }
 187     else{
 188       SetClusterIndex2(i,-3);
 189     }
 190   }
 191   for (Int_t i=0;i<4;i++) {
 192     fDEDX[i] = 0.;
 193     fSDEDX[i] = 1e10;
 194     fNCDEDX[i] = 0;
 195     fNCDEDXInclThres[i] = 0;
 196   }
 197     fDEDX[4] = 0;
 198   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
 199 }
 200
 201 AliTPCseed::AliTPCseed(Double_t xr, Double_t alpha, const Double_t xx[5],
 202                        const Double_t cc[15], Int_t index):
 203   AliTPCtrack(xr, alpha, xx, cc, index),
 204   fEsd(0x0),
 205   fClusterOwner(kFALSE),
 206   fRow(0),
 207   fSector(-1),
 208   fRelativeSector(-1),
 209   fCurrentSigmaY2(-1),
 210   fCurrentSigmaZ2(-1),
 211   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
 212   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
 213   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
 214   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
 215   fErrorY2(1e10),
 216   fErrorZ2(1e10),
 217   fCurrentCluster(0x0),
 218   fCurrentClusterIndex1(-1),
 219   fInDead(kFALSE),
 220   fIsSeeding(kFALSE),
 221   fNoCluster(0),
 222   fSort(0),
 223   fBSigned(kFALSE),
 224   fSeedType(0),
 225   fSeed1(-1),
 226   fSeed2(-1),
 227   fMAngular(0),
 228   fCircular(0),
 229   fPoolID(-1)
 230 {
 231   //
 232   // Constructor
 233   //
 234   fFirstPoint =0;
 235   for (Int_t i=0;i<160;i++) SetClusterIndex2(i,-3);
 236   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
 237   for (Int_t i=0;i<5;i++)   fTPCr[i]=0.2;
 238   for (Int_t i=0;i<4;i++) {
 239     fDEDX[i] = 0.;
 240     fSDEDX[i] = 1e10;
 241     fNCDEDX[i] = 0;
 242     fNCDEDXInclThres[i] = 0;
 243   }
 244     fDEDX[4] = 0;
 245   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
 246 }
 247
 248 AliTPCseed::~AliTPCseed(){
 249   //
 250   // destructor
 251   fNoCluster =0;
 252   if (fClusterOwner){
 253     for (Int_t icluster=0; icluster<160; icluster++){
 254       delete fClusterPointer[icluster];
 255     }
 256   }
 257
 258 }
 259 //_________________________________________________
 260 AliTPCseed & AliTPCseed::operator=(const AliTPCseed &param)
 261 {
 262   //
 263   // assignment operator
 264   // don't touch pool ID
 265   //
 266   if(this!=&param){
 267     AliTPCtrack::operator=(param);
 268     fEsd =param.fEsd;
 269     fClusterOwner = param.fClusterOwner;
 270     if (!fClusterOwner) for(Int_t i = 0;i<160;++i)fClusterPointer[i] = param.fClusterPointer[i];
 271     else                for(Int_t i = 0;i<160;++i) {
 272         delete fClusterPointer[i];
 273         fClusterPointer[i] = new AliTPCclusterMI(*(param.fClusterPointer[i]));
 274       }
 275     // leave out fPoint, they are also not copied in the copy ctor...
 276     // but deleted in the dtor... strange...
 277     fRow            = param.fRow;
 278     fSector         = param.fSector;
 279     fRelativeSector = param.fRelativeSector;
 280     fCurrentSigmaY2 = param.fCurrentSigmaY2;
 281     fCurrentSigmaZ2 = param.fCurrentSigmaZ2;
 282     fErrorY2        = param.fErrorY2;
 283     fErrorZ2        = param.fErrorZ2;
 284     fCurrentCluster = param.fCurrentCluster; // this is not allocated by AliTPCSeed
 285     fCurrentClusterIndex1 = param.fCurrentClusterIndex1;
 286     fInDead         = param.fInDead;
 287     fIsSeeding      = param.fIsSeeding;
 288     fNoCluster      = param.fNoCluster;
 289     fSort           = param.fSort;
 290     fBSigned        = param.fBSigned;
 291     for(Int_t i = 0;i<4;++i){
 292       fDEDX[i]   = param.fDEDX[i];
 293       fSDEDX[i]  = param.fSDEDX[i];
 294       fNCDEDX[i] = param.fNCDEDX[i];
 295       fNCDEDXInclThres[i] = param.fNCDEDXInclThres[i];
 296     }
 297       fDEDX[4]   = param.fDEDX[4];
 298     for(Int_t i = 0;i<AliPID::kSPECIES;++i)fTPCr[i] = param.fTPCr[i];
 299
 300     fSeedType = param.fSeedType;
 301     fSeed1    = param.fSeed1;
 302     fSeed2    = param.fSeed2;
 303     for(Int_t i = 0;i<12;++i)fOverlapLabels[i] = param.fOverlapLabels[i];
 304     fMAngular = param.fMAngular;
 305     fCircular = param.fCircular;
 306     for(int i = 0;i<160;++i)fTrackPoints[i] =  param.fTrackPoints[i];
 307   }
 308   return (*this);
 309 }
 310 //____________________________________________________
 311 AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
 312 {
 313   //
 314   //
 315   return &fTrackPoints[i];
 316 }
 317
 318
 319
 320 Double_t AliTPCseed::GetDensityFirst(Int_t n)
 321 {
 322   //
 323   //
 324   // return cluster for n rows bellow first point
 325   Int_t nfoundable = 1;
 326   Int_t nfound      = 1;
 327   for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
 328     Int_t index = GetClusterIndex2(i);
 329     if (index!=-1) nfoundable++;
 330     if (index>0) nfound++;
 331   }
 332   if (nfoundable<n) return 0;
 333   return Double_t(nfound)/Double_t(nfoundable);
 334
 335 }
 336
 337
 338 void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
 339 {
 340   // get cluster stat.  on given region
 341   //
 342   found       = 0;
 343   foundable   = 0;
 344   shared      =0;
 345   for (Int_t i=first;i<last; i++){
 346     Int_t index = GetClusterIndex2(i);
 347     if (index!=-1) foundable++;
 348     if (index&0x8000) continue;
 349     if (fClusterPointer[i]) {
 350       found++;
 351     }
 352     else
 353       continue;
 354
 355     if (fClusterPointer[i]->IsUsed(10)) {
 356       shared++;
 357       continue;
 358     }
 359     if (!plus2) continue; //take also neighborhoud
 360     //
 361     if ( (i>0) && fClusterPointer[i-1]){
 362       if (fClusterPointer[i-1]->IsUsed(10)) {
 363         shared++;
 364         continue;
 365       }
 366     }
 367     if ( fClusterPointer[i+1]){
 368       if (fClusterPointer[i+1]->IsUsed(10)) {
 369         shared++;
 370         continue;
 371       }
 372     }
 373
 374   }
 375   //if (shared>found){
 376     //Error("AliTPCseed::GetClusterStatistic","problem\n");
 377   //}
 378 }
 379
 380
 381
 382
 383
 384 void AliTPCseed::Reset(Bool_t all)
 385 {
 386   //
 387   //
 388   SetNumberOfClusters(0);
 389   fNFoundable = 0;
 390   SetChi2(0);
 391   ResetCovariance(10.);
 392   /*
 393   if (fTrackPoints){
 394     for (Int_t i=0;i<8;i++){
 395       delete [] fTrackPoints[i];
 396     }
 397     delete fTrackPoints;
 398     fTrackPoints =0;
 399   }
 400   */
 401
 402   if (all){
 403     for (Int_t i=200;i--;) SetClusterIndex2(i,-3);
 404     if (!fClusterOwner) for (Int_t i=160;i--;) fClusterPointer[i]=0;
 405     else                for (Int_t i=160;i--;) {delete fClusterPointer[i]; fClusterPointer[i]=0;}
 406   }
 407
 408 }
 409
 410
 411 void AliTPCseed::Modify(Double_t factor)
 412 {
 413
 414   //------------------------------------------------------------------
 415   //This function makes a track forget its history :)
 416   //------------------------------------------------------------------
 417   if (factor<=0) {
 418     ResetCovariance(10.);
 419     return;
 420   }
 421   ResetCovariance(factor);
 422
 423   SetNumberOfClusters(0);
 424   fNFoundable =0;
 425   SetChi2(0);
 426   fRemoval = 0;
 427   fCurrentSigmaY2 = 0.000005;
 428   fCurrentSigmaZ2 = 0.000005;
 429   fNoCluster     = 0;
 430   //fFirstPoint = 160;
 431   //fLastPoint  = 0;
 432 }
 433
 434
 435
 436
 437 Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
 438 {
 439   //-----------------------------------------------------------------
 440   // This function find proloncation of a track to a reference plane x=xk.
 441   // doesn't change internal state of the track
 442   //-----------------------------------------------------------------
 443
 444   Double_t x1=GetX(), x2=x1+(xk-x1), dx=x2-x1;
 445
 446   if (TMath::Abs(GetSnp()+GetC()*dx) >= AliTPCReconstructor::GetMaxSnpTrack()) {
 447     return 0;
 448   }
 449
 450   //  Double_t y1=fP0, z1=fP1;
 451   Double_t c1=GetSnp(), r1=sqrt((1.-c1)*(1.+c1));
 452   Double_t c2=c1 + GetC()*dx, r2=sqrt((1.-c2)*(1.+c2));
 453
 454   y = GetY();
 455   z = GetZ();
 456   //y += dx*(c1+c2)/(r1+r2);
 457   //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
 458
 459   Double_t dy = dx*(c1+c2)/(r1+r2);
 460   Double_t dz = 0;
 461   //
 462   Double_t delta = GetC()*dx*(c1+c2)/(c1*r2 + c2*r1);
 463   /*
 464   if (TMath::Abs(delta)>0.0001){
 465     dz = fP3*TMath::ASin(delta)/fP4;
 466   }else{
 467     dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
 468   }
 469   */
 470   //  dz =  fP3*AliTPCFastMath::FastAsin(delta)/fP4;
 471   dz =  GetTgl()*TMath::ASin(delta)/GetC();
 472   //
 473   y+=dy;
 474   z+=dz;
 475
 476
 477   return 1;
 478 }
 479
 480
 481 //_____________________________________________________________________________
 482 Double_t AliTPCseed::GetPredictedChi2(const AliCluster *c) const
 483 {
 484   //-----------------------------------------------------------------
 485   // This function calculates a predicted chi2 increment.
 486   //-----------------------------------------------------------------
 487   Double_t p[2]={c->GetY(), c->GetZ()};
 488   Double_t cov[3]={fErrorY2, 0., fErrorZ2};
 489
 490   Float_t dx = ((AliTPCclusterMI*)c)->GetX()-GetX();
 491   if (TMath::Abs(dx)>0){
 492     Float_t ty = TMath::Tan(TMath::ASin(GetSnp()));
 493     Float_t dy = dx*ty;
 494     Float_t dz = dx*TMath::Sqrt(1.+ty*ty)*GetTgl();
 495     p[0] = c->GetY()-dy;
 496     p[1] = c->GetZ()-dz;
 497   }
 498   return AliExternalTrackParam::GetPredictedChi2(p,cov);
 499 }
 500
 501 //_________________________________________________________________________________________
 502
 503
 504 Int_t AliTPCseed::Compare(const TObject *o) const {
 505   //-----------------------------------------------------------------
 506   // This function compares tracks according to the sector - for given sector according z
 507   //-----------------------------------------------------------------
 508   AliTPCseed *t=(AliTPCseed*)o;
 509
 510   if (fSort == 0){
 511     if (t->fRelativeSector>fRelativeSector) return -1;
 512     if (t->fRelativeSector<fRelativeSector) return 1;
 513     Double_t z2 = t->GetZ();
 514     Double_t z1 = GetZ();
 515     if (z2>z1) return 1;
 516     if (z2<z1) return -1;
 517     return 0;
 518   }
 519   else {
 520     Float_t f2 =1;
 521     f2 = 1-20*TMath::Sqrt(t->GetSigma1Pt2())/(t->OneOverPt()+0.0066);
 522     if (t->fBConstrain) f2=1.2;
 523
 524     Float_t f1 =1;
 525     f1 = 1-20*TMath::Sqrt(GetSigma1Pt2())/(OneOverPt()+0.0066);
 526
 527     if (fBConstrain)   f1=1.2;
 528
 529     if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
 530     else return +1;
 531   }
 532 }
 533
 534
 535
 536
 537 //_____________________________________________________________________________
 538 Bool_t AliTPCseed::Update(const AliCluster *c, Double_t chisq, Int_t index)
 539 {
 540   //-----------------------------------------------------------------
 541   // This function associates a cluster with this track.
 542   //-----------------------------------------------------------------
 543   Int_t n=GetNumberOfClusters();
 544   Int_t idx=GetClusterIndex(n);    // save the current cluster index
 545
 546   AliTPCclusterMI cl(*(AliTPCclusterMI*)c);  cl.SetSigmaY2(fErrorY2); cl.SetSigmaZ2(fErrorZ2);
 547
 548   AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
 549
 550   Float_t ty = TMath::Tan(TMath::ASin(GetSnp()));
 551
 552   if(  parcl ){
 553     Int_t padSize = 0;                          // short pads
 554     if (cl.GetDetector() >= 36) {
 555       padSize = 1;                              // medium pads
 556       if (cl.GetRow() > 63) padSize = 2; // long pads
 557     }
 558     Float_t waveCorr = parcl->GetWaveCorrection( padSize, cl.GetZ(), cl.GetMax(),cl.GetPad(), ty );
 559     cl.SetY( cl.GetY() - waveCorr );
 560   }
 561
 562   Float_t dx = ((AliTPCclusterMI*)c)->GetX()-GetX();
 563   if (TMath::Abs(dx)>0){
 564     Float_t dy = dx*ty;
 565     Float_t dz = dx*TMath::Sqrt(1.+ty*ty)*GetTgl();
 566     cl.SetY(cl.GetY()-dy);
 567     cl.SetZ(cl.GetZ()-dz);
 568   }
 569
 570
 571   if (!AliTPCtrack::Update(&cl,chisq,index)) return kFALSE;
 572
 573   if (fCMeanSigmaY2p30<0){
 574     fCMeanSigmaY2p30= c->GetSigmaY2();   //! current mean sigma Y2 - mean30%
 575     fCMeanSigmaZ2p30= c->GetSigmaZ2();   //! current mean sigma Z2 - mean30%
 576     fCMeanSigmaY2p30R = 1;   //! current mean sigma Y2 - mean5%
 577     fCMeanSigmaZ2p30R = 1;   //! current mean sigma Z2 - mean5%
 578   }
 579   //
 580   fCMeanSigmaY2p30= 0.70*fCMeanSigmaY2p30 +0.30*c->GetSigmaY2();
 581   fCMeanSigmaZ2p30= 0.70*fCMeanSigmaZ2p30 +0.30*c->GetSigmaZ2();
 582   if (fCurrentSigmaY2>0){
 583     fCMeanSigmaY2p30R = 0.7*fCMeanSigmaY2p30R  +0.3*c->GetSigmaY2()/fCurrentSigmaY2;
 584     fCMeanSigmaZ2p30R = 0.7*fCMeanSigmaZ2p30R  +0.3*c->GetSigmaZ2()/fCurrentSigmaZ2;
 585   }
 586
 587
 588   SetClusterIndex(n,idx);          // restore the current cluster index
 589   return kTRUE;
 590 }
 591
 592
 593
 594 //_____________________________________________________________________________
 595 Float_t AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t /* onlyused */) {
 596   //-----------------------------------------------------------------
 597   // This funtion calculates dE/dX within the "low" and "up" cuts.
 598   //-----------------------------------------------------------------
 599   // CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal)
 600   AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
 601
 602   Int_t row0 = param->GetNRowLow();
 603   Int_t row1 = row0+param->GetNRowUp1();
 604   Int_t row2 = row1+param->GetNRowUp2();
 605   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
 606   Int_t useTot = 0;
 607   if (recoParam) useTot = (recoParam->GetUseTotCharge())? 0:1;
 608   //
 609   //
 610   TVectorF i1i2;
 611   TVectorF  iro;
 612   TVectorF oro1;
 613   TVectorF oro2;
 614   TVectorF foro;
 615
 616   CookdEdxAnalytical(low,up,useTot ,i1  ,i2,   0, 2, 0, &i1i2);
 617   CookdEdxAnalytical(low,up,useTot ,0   ,row0, 0, 2, 0, &iro);
 618   CookdEdxAnalytical(low,up,useTot ,row0,row1, 0, 2, 0, &oro1);
 619   CookdEdxAnalytical(low,up,useTot ,row1,row2, 0, 2, 0, &oro2);
 620   CookdEdxAnalytical(low,up,useTot ,row0,row2, 0, 2, 0, &foro); // full OROC truncated mean
 621
 622   fDEDX[0]      = i1i2(0);
 623   fDEDX[1]      =  iro(0);
 624   fDEDX[2]      = oro1(0);
 625   fDEDX[3]      = oro2(0);
 626   fDEDX[4]      = foro(0); // full OROC truncated mean
 627   //
 628   fSDEDX[0]     = i1i2(1);
 629   fSDEDX[1]     =  iro(1);
 630   fSDEDX[2]     = oro1(1);
 631   fSDEDX[3]     = oro2(1);
 632   //
 633   fNCDEDX[0]    = TMath::Nint(i1i2(2));
 634   fNCDEDX[1]    = TMath::Nint( iro(2));
 635   fNCDEDX[2]    = TMath::Nint(oro1(2));
 636   fNCDEDX[3]    = TMath::Nint(oro2(2));
 637   //
 638   fNCDEDXInclThres[0]    = TMath::Nint(i1i2(2)+i1i2(9));
 639   fNCDEDXInclThres[1]    = TMath::Nint( iro(2)+ iro(9));
 640   fNCDEDXInclThres[2]    = TMath::Nint(oro1(2)+oro1(9));
 641   fNCDEDXInclThres[3]    = TMath::Nint(oro2(2)+oro2(9));
 642   //
 643   SetdEdx(fDEDX[0]);
 644   return fDEDX[0];
 645
 646 //  return CookdEdxNorm(low,up,0,i1,i2,1,0,2);
 647
 648
 649 //   Float_t amp[200];
 650 //   Float_t angular[200];
 651 //   Float_t weight[200];
 652 //   Int_t index[200];
 653 //   //Int_t nc = 0;
 654 //   Float_t meanlog = 100.;
 655
 656 //   Float_t mean[4]  = {0,0,0,0};
 657 //   Float_t sigma[4] = {1000,1000,1000,1000};
 658 //   Int_t nc[4]      = {0,0,0,0};
 659 //   Float_t norm[4]    = {1000,1000,1000,1000};
 660 //   //
 661 //   //
 662 //   fNShared =0;
 663
 664 //   Float_t gainGG = 1;
 665 //   if (AliTPCcalibDB::Instance()->GetParameters()){
 666 //     gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000.;  //relative gas gain
 667 //   }
 668
 669
 670 //   for (Int_t of =0; of<4; of++){
 671 //     for (Int_t i=of+i1;i<i2;i+=4)
 672 //       {
 673 //      Int_t clindex = fIndex[i];
 674 //      if (clindex<0||clindex&0x8000) continue;
 675
 676 //      //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
 677 //      AliTPCTrackerPoint * point = GetTrackPoint(i);
 678 //      //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
 679 //      //AliTPCTrackerPoint * pointp = 0;
 680 //      //if (i<159) pointp = GetTrackPoint(i+1);
 681
 682 //      if (point==0) continue;
 683 //      AliTPCclusterMI * cl = fClusterPointer[i];
 684 //      if (cl==0) continue;
 685 //      if (onlyused && (!cl->IsUsed(10))) continue;
 686 //      if (cl->IsUsed(11)) {
 687 //        fNShared++;
 688 //        continue;
 689 //      }
 690 //      Int_t   type   = cl->GetType();
 691 //      //if (point->fIsShared){
 692 //      //  fNShared++;
 693 //      //  continue;
 694 //      //}
 695 //      //if (pointm)
 696 //      //  if (pointm->fIsShared) continue;
 697 //      //if (pointp)
 698 //      //  if (pointp->fIsShared) continue;
 699
 700 //      if (type<0) continue;
 701 //      //if (type>10) continue;
 702 //      //if (point->GetErrY()==0) continue;
 703 //      //if (point->GetErrZ()==0) continue;
 704
 705 //      //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
 706 //      //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
 707 //      //if ((ddy*ddy+ddz*ddz)>10) continue;
 708
 709
 710 //      //      if (point->GetCPoint().GetMax()<5) continue;
 711 //      if (cl->GetMax()<5) continue;
 712 //      Float_t angley = point->GetAngleY();
 713 //      Float_t anglez = point->GetAngleZ();
 714
 715 //      Float_t rsigmay2 =  point->GetSigmaY();
 716 //      Float_t rsigmaz2 =  point->GetSigmaZ();
 717 //      /*
 718 //      Float_t ns = 1.;
 719 //      if (pointm){
 720 //        rsigmay +=  pointm->GetTPoint().GetSigmaY();
 721 //        rsigmaz +=  pointm->GetTPoint().GetSigmaZ();
 722 //        ns+=1.;
 723 //      }
 724 //      if (pointp){
 725 //        rsigmay +=  pointp->GetTPoint().GetSigmaY();
 726 //        rsigmaz +=  pointp->GetTPoint().GetSigmaZ();
 727 //        ns+=1.;
 728 //      }
 729 //      rsigmay/=ns;
 730 //      rsigmaz/=ns;
 731 //      */
 732
 733 //      Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
 734
 735 //      Float_t ampc   = 0;     // normalization to the number of electrons
 736 //      if (i>64){
 737 //        //      ampc = 1.*point->GetCPoint().GetMax();
 738 //        ampc = 1.*cl->GetMax();
 739 //        //ampc = 1.*point->GetCPoint().GetQ();
 740 //        //      AliTPCClusterPoint & p = point->GetCPoint();
 741 //        //      Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
 742 //        // Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
 743 //        //Float_t dz =
 744 //        //  TMath::Abs( Int_t(iz) - iz + 0.5);
 745 //        //ampc *= 1.15*(1-0.3*dy);
 746 //        //ampc *= 1.15*(1-0.3*dz);
 747 //        //      Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
 748 //        //ampc               *=zfactor;
 749 //      }
 750 //      else{
 751 //        //ampc = 1.0*point->GetCPoint().GetMax();
 752 //        ampc = 1.0*cl->GetMax();
 753 //        //ampc = 1.0*point->GetCPoint().GetQ();
 754 //        //AliTPCClusterPoint & p = point->GetCPoint();
 755 //        // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
 756 //        //Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
 757 //        //Float_t dz =
 758 //        //  TMath::Abs( Int_t(iz) - iz + 0.5);
 759
 760 //        //ampc *= 1.15*(1-0.3*dy);
 761 //        //ampc *= 1.15*(1-0.3*dz);
 762 //        //    Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
 763 //        //ampc               *=zfactor;
 764
 765 //      }
 766 //      ampc *= 2.0;     // put mean value to channel 50
 767 //      //ampc *= 0.58;     // put mean value to channel 50
 768 //      Float_t w      =  1.;
 769 //      //      if (type>0)  w =  1./(type/2.-0.5);
 770 //      //      Float_t z = TMath::Abs(cl->GetZ());
 771 //      if (i<64) {
 772 //        ampc /= 0.6;
 773 //        //ampc /= (1+0.0008*z);
 774 //      } else
 775 //        if (i>128){
 776 //          ampc /=1.5;
 777 //          //ampc /= (1+0.0008*z);
 778 //        }else{
 779 //          //ampc /= (1+0.0008*z);
 780 //        }
 781
 782 //      if (type<0) {  //amp at the border - lower weight
 783 //        // w*= 2.;
 784
 785 //        continue;
 786 //      }
 787 //      if (rsigma>1.5) ampc/=1.3;  // if big backround
 788 //      amp[nc[of]]        = ampc;
 789 //      amp[nc[of]]       /=gainGG;
 790 //      angular[nc[of]]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
 791 //      weight[nc[of]]     = w;
 792 //      nc[of]++;
 793 //       }
 794
 795 //     TMath::Sort(nc[of],amp,index,kFALSE);
 796 //     Float_t sumamp=0;
 797 //     Float_t sumamp2=0;
 798 //     Float_t sumw=0;
 799 //     //meanlog = amp[index[Int_t(nc[of]*0.33)]];
 800 //     meanlog = 50;
 801 //     for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
 802 //       Float_t ampl      = amp[index[i]]/angular[index[i]];
 803 //       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
 804 //       //
 805 //       sumw    += weight[index[i]];
 806 //       sumamp  += weight[index[i]]*ampl;
 807 //       sumamp2 += weight[index[i]]*ampl*ampl;
 808 //       norm[of]    += angular[index[i]]*weight[index[i]];
 809 //     }
 810 //     if (sumw<1){
 811 //       SetdEdx(0);
 812 //     }
 813 //     else {
 814 //       norm[of] /= sumw;
 815 //       mean[of]  = sumamp/sumw;
 816 //       sigma[of] = sumamp2/sumw-mean[of]*mean[of];
 817 //       if (sigma[of]>0.1)
 818 //      sigma[of] = TMath::Sqrt(sigma[of]);
 819 //       else
 820 //      sigma[of] = 1000;
 821
 822 //     mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
 823 //     //mean  *=(1-0.02*(sigma/(mean*0.17)-1.));
 824 //     //mean *=(1-0.1*(norm-1.));
 825 //     }
 826 //   }
 827
 828 //   Float_t dedx =0;
 829 //   fSdEdx =0;
 830 //   fMAngular =0;
 831 //   //  mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
 832 //   //  mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
 833
 834
 835 //   //  dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
 836 //   //  (  TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
 837
 838 //   Int_t norm2 = 0;
 839 //   Int_t norm3 = 0;
 840 //   for (Int_t i =0;i<4;i++){
 841 //     if (nc[i]>2&&nc[i]<1000){
 842 //       dedx      += mean[i] *nc[i];
 843 //       fSdEdx    += sigma[i]*(nc[i]-2);
 844 //       fMAngular += norm[i] *nc[i];
 845 //       norm2     += nc[i];
 846 //       norm3     += nc[i]-2;
 847 //     }
 848 //     fDEDX[i]  = mean[i];
 849 //     fSDEDX[i] = sigma[i];
 850 //     fNCDEDX[i]= nc[i];
 851 //   }
 852
 853 //   if (norm3>0){
 854 //     dedx   /=norm2;
 855 //     fSdEdx /=norm3;
 856 //     fMAngular/=norm2;
 857 //   }
 858 //   else{
 859 //     SetdEdx(0);
 860 //     return 0;
 861 //   }
 862 //   //  Float_t dedx1 =dedx;
 863 //   /*
 864 //   dedx =0;
 865 //   for (Int_t i =0;i<4;i++){
 866 //     if (nc[i]>2&&nc[i]<1000){
 867 //       mean[i]   = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
 868 //       dedx      += mean[i] *nc[i];
 869 //     }
 870 //     fDEDX[i]  = mean[i];
 871 //   }
 872 //   dedx /= norm2;
 873 //   */
 874
 875
 876 //   SetdEdx(dedx);
 877 //   return dedx;
 878 }
 879
 880 void AliTPCseed::CookPID()
 881 {
 882   //
 883   // cook PID information according dEdx
 884   //
 885   Double_t fRange = 10.;
 886   Double_t fRes   = 0.1;
 887   Double_t fMIP   = 47.;
 888   //
 889   Int_t ns=AliPID::kSPECIES;
 890   Double_t sumr =0;
 891   for (Int_t j=0; j<ns; j++) {
 892     Double_t mass=AliPID::ParticleMass(j);
 893     Double_t mom=GetP();
 894     Double_t dedx=fdEdx/fMIP;
 895     Double_t bethe=AliMathBase::BetheBlochAleph(mom/mass);
 896     Double_t sigma=fRes*bethe;
 897     if (sigma>0.001){
 898       if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 899         fTPCr[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 900         sumr+=fTPCr[j];
 901         continue;
 902       }
 903       fTPCr[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 904       sumr+=fTPCr[j];
 905     }
 906     else{
 907       fTPCr[j]=1.;
 908       sumr+=fTPCr[j];
 909     }
 910   }
 911   for (Int_t j=0; j<ns; j++) {
 912     fTPCr[j]/=sumr;           //normalize
 913   }
 914 }
 915
 916 Double_t AliTPCseed::GetYat(Double_t xk) const {
 917 //-----------------------------------------------------------------
 918 // This function calculates the Y-coordinate of a track at the plane x=xk.
 919 //-----------------------------------------------------------------
 920   if (TMath::Abs(GetSnp())>AliTPCReconstructor::GetMaxSnpTrack()) return 0.; //patch 01 jan 06
 921     Double_t c1=GetSnp(), r1=TMath::Sqrt((1.-c1)*(1.+c1));
 922     Double_t c2=c1+GetC()*(xk-GetX());
 923     if (TMath::Abs(c2)>AliTPCReconstructor::GetMaxSnpTrack()) return 0;
 924     Double_t r2=TMath::Sqrt((1.-c2)*(1.+c2));
 925     return GetY() + (xk-GetX())*(c1+c2)/(r1+r2);
 926 }
 927
 928
 929
 930 Float_t  AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Bool_t shapeNorm,Int_t posNorm, Int_t padNorm, Int_t returnVal){
 931
 932   //
 933   // calculates dedx using the cluster
 934   // low    -  up specify trunc mean range  - default form 0-0.7
 935   // type   -  1 - max charge  or 0- total charge in cluster
 936   //           //2- max no corr 3- total+ correction
 937   // i1-i2  -  the pad-row range used for calculation
 938   // shapeNorm - kTRUE  -taken from OCDB
 939   //
 940   // posNorm   - usage of pos normalization
 941   // padNorm   - pad type normalization
 942   // returnVal - 0 return mean
 943   //           - 1 return RMS
 944   //           - 2 return number of clusters
 945   //
 946   // normalization parametrization taken from AliTPCClusterParam
 947   //
 948   AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
 949   AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
 950   if (!parcl)  return 0;
 951   if (!param) return 0;
 952   Int_t row0 = param->GetNRowLow();
 953   Int_t row1 = row0+param->GetNRowUp1();
 954
 955   Float_t amp[160];
 956   Int_t   indexes[160];
 957   Int_t   ncl=0;
 958   //
 959   //
 960   Float_t gainGG      = 1;  // gas gain factor -always enabled
 961   Float_t gainPad     = 1;  // gain map  - used always
 962   Float_t corrShape   = 1;  // correction due angular effect, diffusion and electron attachment
 963   Float_t corrPos     = 1;  // local position correction - if posNorm enabled
 964   Float_t corrPadType = 1;  // pad type correction - if padNorm enabled
 965   Float_t corrNorm    = 1;  // normalization factor - set Q to channel 50
 966   //
 967   //
 968   //
 969   if (AliTPCcalibDB::Instance()->GetParameters()){
 970     gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000;  //relative gas gain
 971   }
 972
 973   const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
 974   const Float_t kedgey =3.;
 975   //
 976   //
 977   for (Int_t irow=i1; irow<i2; irow++){
 978     AliTPCclusterMI* cluster = GetClusterPointer(irow);
 979     if (!cluster) continue;
 980     if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
 981     Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ();
 982     Int_t  ipad= 0;
 983     if (irow>=row0) ipad=1;
 984     if (irow>=row1) ipad=2;
 985     //
 986     //
 987     //
 988     AliTPCCalPad * gainMap =  AliTPCcalibDB::Instance()->GetDedxGainFactor();
 989     if (gainMap) {
 990       //
 991       // Get gainPad - pad by pad calibration
 992       //
 993       Float_t factor = 1;
 994       AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector());
 995       if (irow < row0) { // IROC
 996         factor = roc->GetValue(irow, TMath::Nint(cluster->GetPad()));
 997       } else {         // OROC
 998         factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad()));
 999       }
1000       if (factor>0.5) gainPad=factor;
1001     }
1002     //
1003     //do position and angular normalization
1004     //
1005     if (shapeNorm){
1006       if (type<=1){
1007         //
1008         AliTPCTrackerPoint * point = GetTrackPoint(irow);
1009         Float_t              ty = TMath::Abs(point->GetAngleY());
1010         Float_t              tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
1011
1012         Float_t dr    = (250.-TMath::Abs(cluster->GetZ()))/250.;
1013         corrShape  = parcl->Qnorm(ipad,type,dr,ty,tz);
1014       }
1015     }
1016
1017     if (posNorm>0){
1018       //
1019       // Do position normalization - relative distance to
1020       // center of pad- time bin
1021       // Work in progress
1022       //      corrPos = parcl->QnormPos(ipad,type, cluster->GetPad(),
1023       //                                cluster->GetTimeBin(), cluster->GetZ(),
1024       //                                cluster->GetSigmaY2(),cluster->GetSigmaZ2(),
1025       //                                cluster->GetMax(),cluster->GetQ());
1026       // scaled response function
1027       Float_t yres0 = parcl->GetRMS0(0,ipad,0,0)/param->GetPadPitchWidth(cluster->GetDetector());
1028       Float_t zres0 = parcl->GetRMS0(1,ipad,0,0)/param->GetZWidth();
1029       //
1030
1031       AliTPCTrackerPoint * point = GetTrackPoint(irow);
1032       Float_t              ty = TMath::Abs(point->GetAngleY());
1033       Float_t              tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
1034
1035       if (type==1) corrPos =
1036         parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
1037                               cluster->GetTimeBin(),ty,tz,yres0,zres0,0.4);
1038       if (type==0) corrPos =
1039         parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
1040                               cluster->GetTimeBin(),ty,tz,yres0,zres0,cluster->GetQ(),2.5,0.4);
1041       if (posNorm==3){
1042         Float_t dr    = (250.-TMath::Abs(cluster->GetZ()))/250.;
1043         Double_t signtgl = (cluster->GetZ()*point->GetAngleZ()>0)? 1:-1;
1044         Double_t p2 = TMath::Abs(TMath::Sin(TMath::ATan(ty)));
1045         Float_t corrHis = parcl->QnormHis(ipad,type,dr,p2,TMath::Abs(point->GetAngleZ())*signtgl);
1046         if (corrHis>0) corrPos*=corrHis;
1047       }
1048
1049     }
1050
1051     if (padNorm==1){
1052       //taken from OCDB
1053       if (type==0 && parcl->QpadTnorm()) corrPadType = (*parcl->QpadTnorm())[ipad];
1054       if (type==1 && parcl->QpadMnorm()) corrPadType = (*parcl->QpadMnorm())[ipad];
1055
1056     }
1057     if (padNorm==2){
1058       corrPadType  =param->GetPadPitchLength(cluster->GetDetector(),cluster->GetRow());
1059       //use hardwired - temp fix
1060       if (type==0) corrNorm=3.;
1061       if (type==1) corrNorm=1.;
1062     }
1063     //
1064     amp[ncl]=charge;
1065     amp[ncl]/=gainGG;
1066     amp[ncl]/=gainPad;
1067     amp[ncl]/=corrShape;
1068     amp[ncl]/=corrPadType;
1069     amp[ncl]/=corrPos;
1070     amp[ncl]/=corrNorm;
1071     //
1072     ncl++;
1073   }
1074
1075   if (type>3) return ncl;
1076   TMath::Sort(ncl,amp, indexes, kFALSE);
1077
1078   if (ncl<10) return 0;
1079
1080   Float_t suma=0;
1081   Float_t suma2=0;
1082   Float_t sumn=0;
1083   Int_t icl0=TMath::Nint(ncl*low);
1084   Int_t icl1=TMath::Nint(ncl*up);
1085   for (Int_t icl=icl0; icl<icl1;icl++){
1086     suma+=amp[indexes[icl]];
1087     suma2+=amp[indexes[icl]]*amp[indexes[icl]];
1088     sumn++;
1089   }
1090   Float_t mean =suma/sumn;
1091   Float_t rms  =TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
1092   //
1093   // do time-dependent correction for pressure and temperature variations
1094   UInt_t runNumber = 1;
1095   Float_t corrTimeGain = 1;
1096   AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
1097   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
1098   if (trans && recoParam->GetUseGainCorrectionTime()>0) {
1099     runNumber = trans->GetCurrentRunNumber();
1100     //AliTPCcalibDB::Instance()->SetRun(runNumber);
1101     TObjArray * timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
1102     if (timeGainSplines) {
1103       UInt_t time = trans->GetCurrentTimeStamp();
1104       AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
1105       AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
1106       if (fitMIP) {
1107         corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time);/*fitMIP->Eval(time);*/
1108       } else {
1109         if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time);/*fitFPcosmic->Eval(time);*/
1110       }
1111     }
1112   }
1113   mean /= corrTimeGain;
1114   rms /= corrTimeGain;
1115   //
1116   if (returnVal==1) return rms;
1117   if (returnVal==2) return ncl;
1118   return mean;
1119 }
1120
1121 Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal, Int_t rowThres, Int_t mode, TVectorT<float> *returnVec){
1122
1123   //
1124   // calculates dedx using the cluster
1125   // low    -  up specify trunc mean range  - default form 0-0.7
1126   // type   -  1 - max charge  or 0- total charge in cluster
1127   //           //2- max no corr 3- total+ correction
1128   // i1-i2  -  the pad-row range used for calculation
1129   //
1130   // posNorm   - usage of pos normalization
1131   // returnVal - 0  return mean
1132   //           - 1  return RMS
1133   //           - 2  return number of clusters
1134   //           - 3  ratio
1135   //           - 4  mean upper half
1136   //           - 5  mean  - lower half
1137   //           - 6  third moment
1138   // mode      - 0 - linear
1139   //           - 1 - logatithmic
1140   // rowThres  - number of rows before and after given pad row to check for clusters below threshold
1141   //
1142   // normalization parametrization taken from AliTPCClusterParam
1143   //
1144   if (returnVec) returnVec->ResizeTo(10);
1145
1146   AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
1147   AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
1148   if (!parcl)  return 0;
1149   if (!param) return 0;
1150   Int_t row0 = param->GetNRowLow();
1151   Int_t row1 = row0+param->GetNRowUp1();
1152
1153   Float_t amp[160];
1154   Int_t   indexes[160];
1155   Int_t   ncl=0;
1156   Int_t   nclBelowThr = 0; // counts number of clusters below threshold
1157   //
1158   //
1159   Float_t gainGG      = 1;  // gas gain factor -always enabled
1160   Float_t gainPad     = 1;  // gain map  - used always
1161   Float_t corrPos     = 1;  // local position correction - if posNorm enabled
1162   //
1163   //
1164   //
1165   if (AliTPCcalibDB::Instance()->GetParameters()){
1166     gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000;  //relative gas gain
1167   }
1168   //
1169   // extract time-dependent correction for pressure and temperature variations
1170   //
1171   UInt_t runNumber = 1;
1172   Float_t corrTimeGain = 1;
1173   TObjArray * timeGainSplines = 0x0;
1174   TGraphErrors * grPadEqual = 0x0;
1175   TGraphErrors*  grChamberGain[3]={0x0,0x0,0x0};
1176   //
1177   AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
1178   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
1179   //
1180   if (recoParam->GetNeighborRowsDedx() == 0) rowThres = 0;
1181   //
1182   if (trans) {
1183       runNumber = trans->GetCurrentRunNumber();
1184       //AliTPCcalibDB::Instance()->SetRun(runNumber);
1185       timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
1186       if (timeGainSplines && recoParam->GetUseGainCorrectionTime()>0) {
1187         UInt_t time = trans->GetCurrentTimeStamp();
1188         AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
1189         AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
1190         if (fitMIP) {
1191           corrTimeGain =  AliTPCcalibDButil::EvalGraphConst(fitMIP, time); /*fitMIP->Eval(time);*/
1192         } else {
1193           if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time); /*fitFPcosmic->Eval(time); */
1194         }
1195         //
1196         if (type==1) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL");
1197         if (type==0) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL");
1198         const char* names[3]={"SHORT","MEDIUM","LONG"};
1199         for (Int_t iPadRegion=0; iPadRegion<3; ++iPadRegion)
1200           grChamberGain[iPadRegion]=(TGraphErrors*)timeGainSplines->FindObject(Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[iPadRegion]));
1201       }
1202   }
1203
1204   const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam
1205   const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
1206   const Float_t kedgey =3.;
1207   //
1208   //
1209   for (Int_t irow=i1; irow<i2; irow++){
1210     AliTPCclusterMI* cluster = GetClusterPointer(irow);
1211     if (!cluster && irow > 1 && irow < 157) {
1212       Bool_t isClBefore = kFALSE;
1213       Bool_t isClAfter  = kFALSE;
1214       for(Int_t ithres = 1; ithres <= rowThres; ithres++) {
1215         AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres);
1216         if (clusterBefore) isClBefore = kTRUE;
1217         AliTPCclusterMI * clusterAfter  = GetClusterPointer(irow + ithres);
1218         if (clusterAfter) isClAfter = kTRUE;
1219       }
1220       if (isClBefore && isClAfter) nclBelowThr++;
1221     }
1222     if (!cluster) continue;
1223     //
1224     //
1225     if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
1226     //
1227     AliTPCTrackerPoint * point = GetTrackPoint(irow);
1228     if (point==0) continue;
1229     Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
1230     if (rsigmay > kClusterShapeCut) continue;
1231     //
1232     if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb
1233     //
1234     Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ();
1235     Int_t  ipad= 0;
1236     if (irow>=row0) ipad=1;
1237     if (irow>=row1) ipad=2;
1238     //
1239     //
1240     //
1241     AliTPCCalPad * gainMap =  AliTPCcalibDB::Instance()->GetDedxGainFactor();
1242     if (gainMap) {
1243       //
1244       // Get gainPad - pad by pad calibration
1245       //
1246       Float_t factor = 1;
1247       AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector());
1248       if (irow < row0) { // IROC
1249         factor = roc->GetValue(irow, TMath::Nint(cluster->GetPad()));
1250       } else {         // OROC
1251         factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad()));
1252       }
1253       if (factor>0.3) gainPad=factor;
1254     }
1255     //
1256     // Do position normalization - relative distance to
1257     // center of pad- time bin
1258
1259     Float_t              ty = TMath::Abs(point->GetAngleY());
1260     Float_t              tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
1261     Float_t yres0 = parcl->GetRMS0(0,ipad,0,0)/param->GetPadPitchWidth(cluster->GetDetector());
1262     Float_t zres0 = parcl->GetRMS0(1,ipad,0,0)/param->GetZWidth();
1263
1264     yres0 *=parcl->GetQnormCorr(ipad, type,0);
1265     zres0 *=parcl->GetQnormCorr(ipad, type,1);
1266     Float_t effLength=parcl->GetQnormCorr(ipad, type,4)*0.5;
1267     Float_t effDiff  =(parcl->GetQnormCorr(ipad, type,2)+parcl->GetQnormCorr(ipad, type,3))*0.5;
1268     //
1269     if (type==1) {
1270       corrPos = parcl->GetQnormCorr(ipad, type,5)*
1271         parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
1272                               cluster->GetTimeBin(),ty,tz,yres0,zres0,effLength,effDiff);
1273       Float_t drm   = 0.5-TMath::Abs(cluster->GetZ()/250.);
1274       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,0)*drm);
1275       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,1)*ty*ty);
1276       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz);
1277       //
1278     }
1279     if (type==0) {
1280       corrPos = parcl->GetQnormCorr(ipad, type,5)*
1281         parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
1282                               cluster->GetTimeBin(),ty,tz,yres0,zres0,cluster->GetQ(),2.5,effLength,effDiff);
1283
1284       Float_t drm   = 0.5-TMath::Abs(cluster->GetZ()/250.);
1285       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,0)*drm);
1286       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,1)*ty*ty);
1287       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz);
1288       //
1289     }
1290     //
1291     // pad region equalization outside of cluster param
1292     //
1293     Float_t gainEqualPadRegion = 1;
1294     if (grPadEqual && recoParam->GetUseGainCorrectionTime()>0) gainEqualPadRegion = grPadEqual->Eval(ipad);
1295     //
1296     // chamber-by-chamber equalization outside gain map
1297     //
1298     Float_t gainChamber = 1;
1299     if (grChamberGain[ipad] && recoParam->GetUseGainCorrectionTime()>0) gainChamber = grChamberGain[ipad]->Eval(cluster->GetDetector());
1300     //
1301     amp[ncl]=charge;
1302     amp[ncl]/=gainGG;
1303     amp[ncl]/=gainPad;
1304     amp[ncl]/=corrPos;
1305     amp[ncl]/=gainEqualPadRegion;
1306     amp[ncl]/=gainChamber;
1307     //
1308     ncl++;
1309   }
1310
1311   if (type==2) return ncl;
1312   TMath::Sort(ncl,amp, indexes, kFALSE);
1313   //
1314   if (ncl<10) return 0;
1315   //
1316   Double_t * ampWithBelow = new Double_t[ncl + nclBelowThr];
1317   for(Int_t iCl = 0; iCl < ncl + nclBelowThr; iCl++) {
1318     if (iCl < nclBelowThr) {
1319       ampWithBelow[iCl] = amp[indexes[0]];
1320     } else {
1321       ampWithBelow[iCl] = amp[indexes[iCl - nclBelowThr]];
1322     }
1323   }
1324   //printf("DEBUG: %i shit %f", nclBelowThr, amp[indexes[0]]);
1325   //
1326   Float_t suma=0;
1327   Float_t suma2=0;
1328   Float_t suma3=0;
1329   Float_t sumaS=0;
1330   Float_t sumn=0;
1331   // upper,and lower part statistic
1332   Float_t sumL=0, sumL2=0, sumLN=0;
1333   Float_t sumD=0, sumD2=0, sumDN=0;
1334
1335   Int_t icl0=TMath::Nint((ncl + nclBelowThr)*low);
1336   Int_t icl1=TMath::Nint((ncl + nclBelowThr)*up);
1337   Int_t iclm=TMath::Nint((ncl + nclBelowThr)*(low +(up+low)*0.5));
1338   //
1339   for (Int_t icl=icl0; icl<icl1;icl++){
1340     if (ampWithBelow[icl]<0.1) continue;
1341     Double_t camp=ampWithBelow[icl]/corrTimeGain;
1342     if (mode==1) camp= TMath::Log(camp);
1343     if (icl<icl1){
1344       suma+=camp;
1345       suma2+=camp*camp;
1346       suma3+=camp*camp*camp;
1347       sumaS+=TMath::Power(TMath::Abs(camp),1./3.);
1348       sumn++;
1349     }
1350     if (icl>iclm){
1351       sumL+=camp;
1352       sumL2+=camp*camp;
1353       sumLN++;
1354       }
1355     if (icl<=iclm){
1356       sumD+=camp;
1357       sumD2+=camp*camp;
1358       sumDN++;
1359     }
1360   }
1361   //
1362   Float_t mean = 0;
1363   Float_t meanL = 0;
1364   Float_t meanD = 0;           // lower half mean
1365   if (sumn > 1e-30)   mean =suma/sumn;
1366   if (sumLN > 1e-30)  meanL =sumL/sumLN;
1367   if (sumDN > 1e-30)  meanD =(sumD/sumDN);
1368   /*
1369   Float_t mean =suma/sumn;
1370   Float_t meanL = sumL/sumLN;
1371   Float_t meanD =(sumD/sumDN);           // lower half mean
1372   */
1373
1374   Float_t rms = 0;
1375   Float_t mean2=0;
1376   Float_t mean3=0;
1377   Float_t meanS=0;
1378
1379   if(sumn>0){
1380     rms = TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
1381     mean2=suma2/sumn;
1382     mean3=suma3/sumn;
1383     meanS=sumaS/sumn;
1384   }
1385
1386   if (mean2>0) mean2=TMath::Power(TMath::Abs(mean2),1./2.);
1387   if (mean3>0) mean3=TMath::Power(TMath::Abs(mean3),1./3.);
1388   if (meanS>0) meanS=TMath::Power(TMath::Abs(meanS),3.);
1389   //
1390   if (mode==1) mean=TMath::Exp(mean);
1391   if (mode==1) meanL=TMath::Exp(meanL);  // upper truncation
1392   if (mode==1) meanD=TMath::Exp(meanD);  // lower truncation
1393   //
1394   delete [] ampWithBelow; //return?
1395
1396
1397   //
1398   if(returnVec){
1399       (*returnVec)(0) = mean;
1400       (*returnVec)(1) = rms;
1401       (*returnVec)(2) = ncl;
1402       (*returnVec)(3) = Double_t(nclBelowThr)/Double_t(nclBelowThr+ncl);
1403       (*returnVec)(4) = meanL;
1404       (*returnVec)(5) = meanD;
1405       (*returnVec)(6) = mean2;
1406       (*returnVec)(7) = mean3;
1407       (*returnVec)(8) = meanS;
1408       (*returnVec)(9) = nclBelowThr;
1409   }
1410
1411   if (returnVal==1) return rms;
1412   if (returnVal==2) return ncl;
1413   if (returnVal==3) return Double_t(nclBelowThr)/Double_t(nclBelowThr+ncl);
1414   if (returnVal==4) return meanL;
1415   if (returnVal==5) return meanD;
1416   if (returnVal==6) return mean2;
1417   if (returnVal==7) return mean3;
1418   if (returnVal==8) return meanS;
1419   if (returnVal==9) return nclBelowThr;
1420   return mean;
1421 }
1422
1423
1424
1425
1426 Float_t  AliTPCseed::CookShape(Int_t type){
1427   //
1428   //
1429   //
1430  //-----------------------------------------------------------------
1431   // This funtion calculates dE/dX within the "low" and "up" cuts.
1432   //-----------------------------------------------------------------
1433   Float_t means=0;
1434   Float_t meanc=0;
1435   for (Int_t i =0; i<160;i++)    {
1436     AliTPCTrackerPoint * point = GetTrackPoint(i);
1437     if (point==0) continue;
1438
1439     AliTPCclusterMI * cl = fClusterPointer[i];
1440     if (cl==0) continue;
1441
1442     Float_t rsigmay =  TMath::Sqrt(point->GetSigmaY());
1443     Float_t rsigmaz =  TMath::Sqrt(point->GetSigmaZ());
1444     Float_t rsigma =   (rsigmay+rsigmaz)*0.5;
1445     if (type==0) means+=rsigma;
1446     if (type==1) means+=rsigmay;
1447     if (type==2) means+=rsigmaz;
1448     meanc++;
1449   }
1450   Float_t mean = (meanc>0)? means/meanc:0;
1451   return mean;
1452 }
1453
1454
1455
1456 Int_t  AliTPCseed::RefitTrack(AliTPCseed *seed, AliExternalTrackParam * parin, AliExternalTrackParam * parout){
1457   //
1458   // Refit the track
1459   // return value - number of used clusters
1460   //
1461   //
1462   const Int_t kMinNcl =10;
1463   AliTPCseed *track=new AliTPCseed(*seed);
1464   Int_t sector=-1;
1465   // reset covariance
1466   //
1467   Double_t covar[15];
1468   for (Int_t i=0;i<15;i++) covar[i]=0;
1469   covar[0]=10.*10.;
1470   covar[2]=10.*10.;
1471   covar[5]=10.*10./(64.*64.);
1472   covar[9]=10.*10./(64.*64.);
1473   covar[14]=1*1;
1474   //
1475
1476   Float_t xmin=1000, xmax=-10000;
1477   Int_t imin=158, imax=0;
1478   for (Int_t i=0;i<160;i++) {
1479     AliTPCclusterMI *c=track->GetClusterPointer(i);
1480     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue;
1481     if (sector<0) sector = c->GetDetector();
1482     if (c->GetX()<xmin) xmin=c->GetX();
1483     if (c->GetX()>xmax) xmax=c->GetX();
1484     if (i<imin) imin=i;
1485     if (i>imax) imax=i;
1486   }
1487   if(imax-imin<kMinNcl) {
1488     delete track;
1489     return 0 ;
1490   }
1491   // Not succes to rotate
1492   if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1493     delete track;
1494     return 0;
1495   }
1496   //
1497   //
1498   // fit from inner to outer row
1499   //
1500   AliExternalTrackParam paramIn;
1501   AliExternalTrackParam paramOut;
1502   Bool_t isOK=kTRUE;
1503   Int_t ncl=0;
1504   //
1505   //
1506   //
1507   for (Int_t i=imin; i<=imax; i++){
1508     AliTPCclusterMI *c=track->GetClusterPointer(i);
1509     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue;
1510     //    if (RejectCluster(c,track)) continue;
1511     sector = (c->GetDetector()%18);
1512     if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1513       //continue;
1514     }
1515     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1516     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1517     if (!track->PropagateTo(r[0])) {
1518       isOK=kFALSE;
1519     }
1520     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1521   }
1522   if (!isOK) { delete track; return 0;}
1523   track->AddCovariance(covar);
1524   //
1525   //
1526   //
1527   for (Int_t i=imax; i>=imin; i--){
1528     AliTPCclusterMI *c=track->GetClusterPointer(i);
1529     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue;
1530     //if (RejectCluster(c,track)) continue;
1531     sector = (c->GetDetector()%18);
1532     if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1533       //continue;
1534     }
1535     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1536     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1537     if (!track->PropagateTo(r[0])) {
1538       isOK=kFALSE;
1539     }
1540     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1541   }
1542   //if (!isOK) { delete track; return 0;}
1543   paramIn = *track;
1544   track->AddCovariance(covar);
1545   //
1546   //
1547   for (Int_t i=imin; i<=imax; i++){
1548     AliTPCclusterMI *c=track->GetClusterPointer(i);
1549     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue;
1550     sector = (c->GetDetector()%18);
1551     if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1552       //continue;
1553     }
1554     ncl++;
1555     //if (RejectCluster(c,track)) continue;
1556     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1557     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1558     if (!track->PropagateTo(r[0])) {
1559       isOK=kFALSE;
1560     }
1561     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1562   }
1563   //if (!isOK) { delete track; return 0;}
1564   paramOut=*track;
1565   //
1566   //
1567   //
1568   if (parin) (*parin)=paramIn;
1569   if (parout) (*parout)=paramOut;
1570   delete track;
1571   return ncl;
1572 }
1573
1574
1575
1576 Bool_t AliTPCseed::RefitTrack(AliTPCseed* /*seed*/, Bool_t /*out*/){
1577   //
1578   //
1579   //
1580   return kFALSE;
1581 }
1582
1583
1584
1585
1586
1587
1588 void  AliTPCseed::GetError(AliTPCclusterMI* cluster, AliExternalTrackParam * param,
1589                                   Double_t& erry, Double_t &errz)
1590 {
1591   //
1592   // Get cluster error at given position
1593   //
1594   AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam();
1595   Double_t tany,tanz;
1596   Double_t snp1=param->GetSnp();
1597   tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1598   //
1599   Double_t tgl1=param->GetTgl();
1600   tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1601   //
1602   Int_t padSize = 0;                          // short pads
1603   if (cluster->GetDetector() >= 36) {
1604     padSize = 1;                              // medium pads
1605     if (cluster->GetRow() > 63) padSize = 2; // long pads
1606   }
1607
1608   erry  = clusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany) );
1609   errz  = clusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) );
1610 }
1611
1612
1613 void  AliTPCseed::GetShape(AliTPCclusterMI* cluster, AliExternalTrackParam * param,
1614                                   Double_t& rmsy, Double_t &rmsz)
1615 {
1616   //
1617   // Get cluster error at given position
1618   //
1619   AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam();
1620   Double_t tany,tanz;
1621   Double_t snp1=param->GetSnp();
1622   tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1623   //
1624   Double_t tgl1=param->GetTgl();
1625   tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1626   //
1627   Int_t padSize = 0;                          // short pads
1628   if (cluster->GetDetector() >= 36) {
1629     padSize = 1;                              // medium pads
1630     if (cluster->GetRow() > 63) padSize = 2; // long pads
1631   }
1632
1633   rmsy  = clusterParam->GetRMSQ( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany), TMath::Abs(cluster->GetMax()) );
1634   rmsz  = clusterParam->GetRMSQ( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) ,TMath::Abs(cluster->GetMax()));
1635 }
1636
1637
1638
1639 Double_t AliTPCseed::GetQCorrGeom(Float_t ty, Float_t tz){
1640   //Geoetrical
1641   //ty    - tangent in local y direction
1642   //tz    - tangent
1643   //
1644   Float_t norm=TMath::Sqrt(1+ty*ty+tz*tz);
1645   return norm;
1646 }
1647
1648 Double_t AliTPCseed::GetQCorrShape(Int_t ipad, Int_t type,Float_t z, Float_t ty, Float_t tz, Float_t /*q*/, Float_t /*thr*/){
1649   //
1650   // Q normalization
1651   //
1652   // return value =  Q Normalization factor
1653   // Normalization - 1 - shape factor part for full drift
1654   //                 1 - electron attachment for 0 drift
1655
1656   // Input parameters:
1657   //
1658   // ipad - 0 short pad
1659   //        1 medium pad
1660   //        2 long pad
1661   //
1662   // type - 0 qmax
1663   //      - 1 qtot
1664   //
1665   //z     - z position (-250,250 cm)
1666   //ty    - tangent in local y direction
1667   //tz    - tangent
1668   //
1669
1670   AliTPCClusterParam * paramCl = AliTPCcalibDB::Instance()->GetClusterParam();
1671   AliTPCParam   * paramTPC = AliTPCcalibDB::Instance()->GetParameters();
1672
1673   if (!paramCl) return 1;
1674   //
1675   Double_t dr =  250.-TMath::Abs(z);
1676   Double_t sy =  paramCl->GetRMS0( 0,ipad, dr, TMath::Abs(ty));
1677   Double_t sy0=  paramCl->GetRMS0(0,ipad, 250, 0);
1678   Double_t sz =  paramCl->GetRMS0( 1,ipad, dr, TMath::Abs(tz));
1679   Double_t sz0=  paramCl->GetRMS0(1,ipad, 250, 0);
1680
1681   Double_t sfactorMax = TMath::Sqrt(sy0*sz0/(sy*sz));
1682
1683
1684   Double_t dt = 1000000*(dr/paramTPC->GetDriftV());  //time in microsecond
1685   Double_t attProb = TMath::Exp(-paramTPC->GetAttCoef()*paramTPC->GetOxyCont()*dt);
1686   //
1687   //
1688   if (type==0) return sfactorMax*attProb;
1689
1690   return attProb;
1691
1692
1693 }
1694
1695
1696 /*
1697 //_______________________________________________________________________
1698 Float_t AliTPCseed::GetTPCClustInfo(Int_t nNeighbours, Int_t type, Int_t row0, Int_t row1, TVectorT<float> *returnVec)
1699 {
1700   //
1701   // TPC cluster information
1702   // type 0: get fraction of found/findable clusters with neighbourhood definition
1703   //      1: found clusters
1704   //      2: findable (number of clusters above and below threshold)
1705   //
1706   // definition of findable clusters:
1707   //            a cluster is defined as findable if there is another cluster
1708   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1709   //           effects with a very simple algorithm.
1710   //
1711
1712   const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam
1713   const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
1714   const Float_t kedgey =3.;
1715
1716   Float_t ncl = 0;
1717   Float_t nclBelowThr = 0; // counts number of clusters below threshold
1718
1719   for (Int_t irow=row0; irow<row1; irow++){
1720     AliTPCclusterMI* cluster = GetClusterPointer(irow);
1721
1722     if (!cluster && irow > 1 && irow < 157) {
1723       Bool_t isClBefore = kFALSE;
1724       Bool_t isClAfter  = kFALSE;
1725       for(Int_t ithres = 1; ithres <= nNeighbours; ithres++) {
1726         AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres);
1727         if (clusterBefore) isClBefore = kTRUE;
1728         AliTPCclusterMI * clusterAfter  = GetClusterPointer(irow + ithres);
1729         if (clusterAfter) isClAfter = kTRUE;
1730       }
1731       if (isClBefore && isClAfter) nclBelowThr++;
1732     }
1733     if (!cluster) continue;
1734     //
1735     //
1736     if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
1737     //
1738     AliTPCTrackerPoint * point = GetTrackPoint(irow);
1739     if (point==0) continue;
1740     Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
1741     if (rsigmay > kClusterShapeCut) continue;
1742     //
1743     if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb
1744     ncl++;
1745   }
1746   if(returnVec->GetNoElements != 3){
1747       returnVec->ResizeTo(3);
1748   }
1749   Float_t nclAll = nclBelowThr+ncl;
1750   returnVec(0) = nclAll>0?ncl/nclAll:0;
1751   returnVec(1) = ncl;
1752   returnVec(2) = nclAll;
1753
1754   if(ncl<10)
1755     return 0;
1756   if(type==0)
1757     if(nclAll>0)
1758       return ncl/nclAll;
1759   if(type==1)
1760     return ncl;
1761   if(type==2)
1762     return nclAll;
1763   return 0;
1764 }
1765 */
1766 //_______________________________________________________________________
1767 Int_t AliTPCseed::GetNumberOfClustersIndices() {
1768   Int_t ncls = 0;
1769   for (int i=0; i < 160; i++) {
1770     if ((fIndex[i] & 0x8000) == 0)
1771       ncls++;
1772   }
1773   return ncls;
1774 }
1775
1776 //_______________________________________________________________________
1777 void AliTPCseed::Clear(Option_t*)
1778 {
1779   // formally seed may allocate memory for clusters (althought this should not happen for
1780   // the seeds in the pool). Hence we need this method for fwd. compatibility
1781   if (fClusterOwner) for (int i=160;i--;) {delete fClusterPointer[i]; fClusterPointer[i] = 0;}
1782 }
1783
1784 TObject* AliTPCseed::Clone(const char* /*newname*/) const
1785 {
1786   // temporary override TObject::Clone to avoid crashes in reco
1787   AliTPCseed* src = (AliTPCseed*)this;
1788   AliTPCseed* dst = new AliTPCseed(*src,fClusterOwner);
1789   return dst;
1790 }