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