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     //
1086     amp[ncl]=charge;
1087     amp[ncl]/=gainGG;                 // normalized gas gain
1088     amp[ncl]/=gainPad;                //
1089     amp[ncl]/=corrShape;
1090     amp[ncl]/=corrPadType;
1091     amp[ncl]/=corrPos;
1092     amp[ncl]/=corrNorm;
1093     //
1094     ncl++;
1095   }
1096
1097   if (type>3) return ncl;
1098   TMath::Sort(ncl,amp, indexes, kFALSE);
1099
1100   if (ncl<10) return 0;
1101
1102   Float_t suma=0;
1103   Float_t suma2=0;
1104   Float_t sumn=0;
1105   Int_t icl0=TMath::Nint(ncl*low);
1106   Int_t icl1=TMath::Nint(ncl*up);
1107   for (Int_t icl=icl0; icl<icl1;icl++){
1108     suma+=amp[indexes[icl]];
1109     suma2+=amp[indexes[icl]]*amp[indexes[icl]];
1110     sumn++;
1111   }
1112   Float_t mean =suma/sumn;
1113   Float_t rms  =TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
1114   //
1115   // do time-dependent correction for pressure and temperature variations
1116   UInt_t runNumber = 1;
1117   Float_t corrTimeGain = 1;
1118   AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
1119   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
1120   if (trans && recoParam->GetUseGainCorrectionTime()>0) {
1121     runNumber = trans->GetCurrentRunNumber();
1122     //AliTPCcalibDB::Instance()->SetRun(runNumber);
1123     TObjArray * timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
1124     if (timeGainSplines) {
1125       UInt_t time = trans->GetCurrentTimeStamp();
1126       AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
1127       AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
1128       if (fitMIP) {
1129         corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time);/*fitMIP->Eval(time);*/
1130       } else {
1131         if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time);/*fitFPcosmic->Eval(time);*/
1132       }
1133     }
1134   }
1135   mean /= corrTimeGain;
1136   rms /= corrTimeGain;
1137   //
1138   if (returnVal==1) return rms;
1139   if (returnVal==2) return ncl;
1140   return mean;
1141 }
1142
1143 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){
1144
1145   //
1146   // calculates dedx using the cluster
1147   // low    -  up specify trunc mean range  - default form 0-0.7
1148   // type   -  1 - max charge  or 0- total charge in cluster
1149   //           //2- max no corr 3- total+ correction
1150   // i1-i2  -  the pad-row range used for calculation
1151   //
1152   // posNorm   - usage of pos normalization
1153   // returnVal - 0  return mean
1154   //           - 1  return RMS
1155   //           - 2  return number of clusters
1156   //           - 3  ratio
1157   //           - 4  mean upper half
1158   //           - 5  mean  - lower half
1159   //           - 6  third moment
1160   // mode      - 0 - linear
1161   //           - 1 - logatithmic
1162   // rowThres  - number of rows before and after given pad row to check for clusters below threshold
1163   //
1164   // normalization parametrization taken from AliTPCClusterParam
1165   //
1166   if (returnVec) returnVec->ResizeTo(10);
1167
1168   AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
1169   AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
1170   if (!parcl)  return 0;
1171   if (!param) return 0;
1172   Int_t row0 = param->GetNRowLow();
1173   Int_t row1 = row0+param->GetNRowUp1();
1174
1175   Float_t amp[160];
1176   Int_t   indexes[160];
1177   Int_t   ncl=0;
1178   Int_t   nclBelowThr = 0; // counts number of clusters below threshold
1179   //
1180   //
1181   Float_t gainGG      = 1;  // gas gain factor -always enabled
1182   Float_t gainPad     = 1;  // gain map  - used always
1183   Float_t corrPos     = 1;  // local position correction - if posNorm enabled
1184   //
1185   //
1186   //
1187   if (AliTPCcalibDB::Instance()->GetParameters()){
1188     gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000;  //relative gas gain
1189   }
1190   //
1191   // extract time-dependent correction for pressure and temperature variations
1192   //
1193   UInt_t runNumber = 1;
1194   Float_t corrTimeGain = 1;
1195   TObjArray * timeGainSplines = 0x0;
1196   TGraphErrors * grPadEqual = 0x0;
1197   TGraphErrors*  grChamberGain[3]={0x0,0x0,0x0};
1198   //
1199   AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
1200   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
1201   //
1202   if (recoParam->GetNeighborRowsDedx() == 0) rowThres = 0;
1203   UInt_t time = 1;//
1204   //
1205   if (trans) {
1206       runNumber = trans->GetCurrentRunNumber();
1207       //AliTPCcalibDB::Instance()->SetRun(runNumber);
1208       timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
1209       if (timeGainSplines && recoParam->GetUseGainCorrectionTime()>0) {
1210         time = trans->GetCurrentTimeStamp();
1211         AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
1212         AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
1213         if (fitMIP) {
1214           corrTimeGain =  AliTPCcalibDButil::EvalGraphConst(fitMIP, time); /*fitMIP->Eval(time);*/
1215         } else {
1216           if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time); /*fitFPcosmic->Eval(time); */
1217         }
1218         //
1219         if (type==1) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL");
1220         if (type==0) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL");
1221         const char* names[3]={"SHORT","MEDIUM","LONG"};
1222         for (Int_t iPadRegion=0; iPadRegion<3; ++iPadRegion)
1223           grChamberGain[iPadRegion]=(TGraphErrors*)timeGainSplines->FindObject(Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[iPadRegion]));
1224       }
1225   }
1226
1227   const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam
1228   const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
1229   const Float_t kedgey =3.;
1230   //
1231   //
1232   for (Int_t irow=i1; irow<i2; irow++){
1233     AliTPCclusterMI* cluster = GetClusterPointer(irow);
1234     if (!cluster && irow > 1 && irow < 157) {
1235       Bool_t isClBefore = kFALSE;
1236       Bool_t isClAfter  = kFALSE;
1237       for(Int_t ithres = 1; ithres <= rowThres; ithres++) {
1238         AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres);
1239         if (clusterBefore) isClBefore = kTRUE;
1240         AliTPCclusterMI * clusterAfter  = GetClusterPointer(irow + ithres);
1241         if (clusterAfter) isClAfter = kTRUE;
1242       }
1243       if (isClBefore && isClAfter) nclBelowThr++;
1244     }
1245     if (!cluster) continue;
1246     //
1247     //
1248     if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
1249     //
1250     AliTPCTrackerPoint * point = GetTrackPoint(irow);
1251     if (point==0) continue;
1252     Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
1253     if (rsigmay > kClusterShapeCut) continue;
1254     //
1255     if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb
1256     //
1257     Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ();
1258     Int_t  ipad= 0;
1259     if (irow>=row0) ipad=1;
1260     if (irow>=row1) ipad=2;
1261     //
1262     //
1263     //
1264     AliTPCCalPad * gainMap =  AliTPCcalibDB::Instance()->GetDedxGainFactor();
1265     if (gainMap) {
1266       //
1267       // Get gainPad - pad by pad calibration
1268       //
1269       Float_t factor = 1;
1270       AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector());
1271       if (irow < row0) { // IROC
1272         factor = roc->GetValue(irow, TMath::Nint(cluster->GetPad()));
1273       } else {         // OROC
1274         factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad()));
1275       }
1276       if (factor>0.3) gainPad=factor;
1277     }
1278     //
1279     // Do position normalization - relative distance to
1280     // center of pad- time bin
1281
1282     Float_t              ty = TMath::Abs(point->GetAngleY());
1283     Float_t              tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
1284     Float_t yres0 = parcl->GetRMS0(0,ipad,0,0)/param->GetPadPitchWidth(cluster->GetDetector());
1285     Float_t zres0 = parcl->GetRMS0(1,ipad,0,0)/param->GetZWidth();
1286
1287     yres0 *=parcl->GetQnormCorr(ipad, type,0);
1288     zres0 *=parcl->GetQnormCorr(ipad, type,1);
1289     Float_t effLength=parcl->GetQnormCorr(ipad, type,4)*0.5;
1290     Float_t effDiff  =(parcl->GetQnormCorr(ipad, type,2)+parcl->GetQnormCorr(ipad, type,3))*0.5;
1291     //
1292     if (type==1) {
1293       corrPos = parcl->GetQnormCorr(ipad, type,5)*
1294         parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
1295                               cluster->GetTimeBin(),ty,tz,yres0,zres0,effLength,effDiff);
1296       Float_t drm   = 0.5-TMath::Abs(cluster->GetZ()/250.);
1297       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,0)*drm);
1298       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,1)*ty*ty);
1299       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz);
1300       //
1301     }
1302     if (type==0) {
1303       corrPos = parcl->GetQnormCorr(ipad, type,5)*
1304         parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
1305                               cluster->GetTimeBin(),ty,tz,yres0,zres0,cluster->GetQ(),2.5,effLength,effDiff);
1306
1307       Float_t drm   = 0.5-TMath::Abs(cluster->GetZ()/250.);
1308       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,0)*drm);
1309       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,1)*ty*ty);
1310       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz);
1311       //
1312     }
1313     //
1314     // pad region equalization outside of cluster param
1315     //
1316     Float_t gainEqualPadRegion = 1;
1317     if (grPadEqual && recoParam->GetUseGainCorrectionTime()>0) gainEqualPadRegion = grPadEqual->Eval(ipad);
1318     //
1319     // chamber-by-chamber equalization outside gain map
1320     //
1321     Float_t gainChamber = 1;
1322     if (grChamberGain[ipad] && recoParam->GetUseGainCorrectionTime()>0) {
1323       gainChamber = grChamberGain[ipad]->Eval(cluster->GetDetector());
1324       if (gainChamber==0) gainChamber=1; // in case old calibation was used before use no correction
1325     }
1326     //
1327     //
1328     //
1329     Double_t correctionHVandPT = AliTPCcalibDB::Instance()->GetGainCorrectionHVandPT(time, runNumber,cluster->GetDetector(), 5 , recoParam->GetGainCorrectionHVandPTMode());
1330     //
1331     amp[ncl]=charge;
1332     amp[ncl]/=gainGG;               // nominal gas gain
1333     amp[ncl]/=correctionHVandPT;    // correction for the HV and P/T - time dependent
1334     amp[ncl]/=gainPad;              //
1335     amp[ncl]/=corrPos;
1336     amp[ncl]/=gainEqualPadRegion;
1337     amp[ncl]/=gainChamber;
1338     //
1339     ncl++;
1340   }
1341
1342   if (type==2) return ncl;
1343   TMath::Sort(ncl,amp, indexes, kFALSE);
1344   //
1345   if (ncl<10) return 0;
1346   //
1347   Double_t * ampWithBelow = new Double_t[ncl + nclBelowThr];
1348   for(Int_t iCl = 0; iCl < ncl + nclBelowThr; iCl++) {
1349     if (iCl < nclBelowThr) {
1350       ampWithBelow[iCl] = amp[indexes[0]];
1351     } else {
1352       ampWithBelow[iCl] = amp[indexes[iCl - nclBelowThr]];
1353     }
1354   }
1355   //printf("DEBUG: %i shit %f", nclBelowThr, amp[indexes[0]]);
1356   //
1357   Float_t suma=0;
1358   Float_t suma2=0;
1359   Float_t suma3=0;
1360   Float_t sumaS=0;
1361   Float_t sumn=0;
1362   // upper,and lower part statistic
1363   Float_t sumL=0, sumL2=0, sumLN=0;
1364   Float_t sumD=0, sumD2=0, sumDN=0;
1365
1366   Int_t icl0=TMath::Nint((ncl + nclBelowThr)*low);
1367   Int_t icl1=TMath::Nint((ncl + nclBelowThr)*up);
1368   Int_t iclm=TMath::Nint((ncl + nclBelowThr)*(low +(up+low)*0.5));
1369   //
1370   for (Int_t icl=icl0; icl<icl1;icl++){
1371     if (ampWithBelow[icl]<0.1) continue;
1372     Double_t camp=ampWithBelow[icl]/corrTimeGain;
1373     if (mode==1) camp= TMath::Log(camp);
1374     if (icl<icl1){
1375       suma+=camp;
1376       suma2+=camp*camp;
1377       suma3+=camp*camp*camp;
1378       sumaS+=TMath::Power(TMath::Abs(camp),1./3.);
1379       sumn++;
1380     }
1381     if (icl>iclm){
1382       sumL+=camp;
1383       sumL2+=camp*camp;
1384       sumLN++;
1385       }
1386     if (icl<=iclm){
1387       sumD+=camp;
1388       sumD2+=camp*camp;
1389       sumDN++;
1390     }
1391   }
1392   //
1393   Float_t mean = 0;
1394   Float_t meanL = 0;
1395   Float_t meanD = 0;           // lower half mean
1396   if (sumn > 1e-30)   mean =suma/sumn;
1397   if (sumLN > 1e-30)  meanL =sumL/sumLN;
1398   if (sumDN > 1e-30)  meanD =(sumD/sumDN);
1399   /*
1400   Float_t mean =suma/sumn;
1401   Float_t meanL = sumL/sumLN;
1402   Float_t meanD =(sumD/sumDN);           // lower half mean
1403   */
1404
1405   Float_t rms = 0;
1406   Float_t mean2=0;
1407   Float_t mean3=0;
1408   Float_t meanS=0;
1409
1410   if(sumn>0){
1411     rms = TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
1412     mean2=suma2/sumn;
1413     mean3=suma3/sumn;
1414     meanS=sumaS/sumn;
1415   }
1416
1417   if (mean2>0) mean2=TMath::Power(TMath::Abs(mean2),1./2.);
1418   if (mean3>0) mean3=TMath::Power(TMath::Abs(mean3),1./3.);
1419   if (meanS>0) meanS=TMath::Power(TMath::Abs(meanS),3.);
1420   //
1421   if (mode==1) mean=TMath::Exp(mean);
1422   if (mode==1) meanL=TMath::Exp(meanL);  // upper truncation
1423   if (mode==1) meanD=TMath::Exp(meanD);  // lower truncation
1424   //
1425   delete [] ampWithBelow; //return?
1426
1427
1428   //
1429   if(returnVec){
1430       (*returnVec)(0) = mean;
1431       (*returnVec)(1) = rms;
1432       (*returnVec)(2) = ncl;
1433       (*returnVec)(3) = Double_t(nclBelowThr)/Double_t(nclBelowThr+ncl);
1434       (*returnVec)(4) = meanL;
1435       (*returnVec)(5) = meanD;
1436       (*returnVec)(6) = mean2;
1437       (*returnVec)(7) = mean3;
1438       (*returnVec)(8) = meanS;
1439       (*returnVec)(9) = nclBelowThr;
1440   }
1441
1442   if (returnVal==1) return rms;
1443   if (returnVal==2) return ncl;
1444   if (returnVal==3) return Double_t(nclBelowThr)/Double_t(nclBelowThr+ncl);
1445   if (returnVal==4) return meanL;
1446   if (returnVal==5) return meanD;
1447   if (returnVal==6) return mean2;
1448   if (returnVal==7) return mean3;
1449   if (returnVal==8) return meanS;
1450   if (returnVal==9) return nclBelowThr;
1451   return mean;
1452 }
1453
1454
1455
1456
1457 Float_t  AliTPCseed::CookShape(Int_t type){
1458   //
1459   //
1460   //
1461  //-----------------------------------------------------------------
1462   // This funtion calculates dE/dX within the "low" and "up" cuts.
1463   //-----------------------------------------------------------------
1464   Float_t means=0;
1465   Float_t meanc=0;
1466   for (Int_t i =0; i<160;i++)    {
1467     AliTPCTrackerPoint * point = GetTrackPoint(i);
1468     if (point==0) continue;
1469
1470     AliTPCclusterMI * cl = fClusterPointer[i];
1471     if (cl==0) continue;
1472
1473     Float_t rsigmay =  TMath::Sqrt(point->GetSigmaY());
1474     Float_t rsigmaz =  TMath::Sqrt(point->GetSigmaZ());
1475     Float_t rsigma =   (rsigmay+rsigmaz)*0.5;
1476     if (type==0) means+=rsigma;
1477     if (type==1) means+=rsigmay;
1478     if (type==2) means+=rsigmaz;
1479     meanc++;
1480   }
1481   Float_t mean = (meanc>0)? means/meanc:0;
1482   return mean;
1483 }
1484
1485
1486
1487 Int_t  AliTPCseed::RefitTrack(AliTPCseed *seed, AliExternalTrackParam * parin, AliExternalTrackParam * parout){
1488   //
1489   // Refit the track
1490   // return value - number of used clusters
1491   //
1492   //
1493   const Int_t kMinNcl =10;
1494   AliTPCseed *track=new AliTPCseed(*seed);
1495   Int_t sector=-1;
1496   // reset covariance
1497   //
1498   Double_t covar[15];
1499   for (Int_t i=0;i<15;i++) covar[i]=0;
1500   covar[0]=10.*10.;
1501   covar[2]=10.*10.;
1502   covar[5]=10.*10./(64.*64.);
1503   covar[9]=10.*10./(64.*64.);
1504   covar[14]=1*1;
1505   //
1506
1507   Float_t xmin=1000, xmax=-10000;
1508   Int_t imin=158, imax=0;
1509   for (Int_t i=0;i<160;i++) {
1510     AliTPCclusterMI *c=track->GetClusterPointer(i);
1511     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue;
1512     if (sector<0) sector = c->GetDetector();
1513     if (c->GetX()<xmin) xmin=c->GetX();
1514     if (c->GetX()>xmax) xmax=c->GetX();
1515     if (i<imin) imin=i;
1516     if (i>imax) imax=i;
1517   }
1518   if(imax-imin<kMinNcl) {
1519     delete track;
1520     return 0 ;
1521   }
1522   // Not succes to rotate
1523   if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1524     delete track;
1525     return 0;
1526   }
1527   //
1528   //
1529   // fit from inner to outer row
1530   //
1531   AliExternalTrackParam paramIn;
1532   AliExternalTrackParam paramOut;
1533   Bool_t isOK=kTRUE;
1534   Int_t ncl=0;
1535   //
1536   //
1537   //
1538   for (Int_t i=imin; i<=imax; i++){
1539     AliTPCclusterMI *c=track->GetClusterPointer(i);
1540     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue;
1541     //    if (RejectCluster(c,track)) continue;
1542     sector = (c->GetDetector()%18);
1543     if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1544       //continue;
1545     }
1546     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1547     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1548     if (!track->PropagateTo(r[0])) {
1549       isOK=kFALSE;
1550     }
1551     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1552   }
1553   if (!isOK) { delete track; return 0;}
1554   track->AddCovariance(covar);
1555   //
1556   //
1557   //
1558   for (Int_t i=imax; i>=imin; i--){
1559     AliTPCclusterMI *c=track->GetClusterPointer(i);
1560     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue;
1561     //if (RejectCluster(c,track)) continue;
1562     sector = (c->GetDetector()%18);
1563     if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1564       //continue;
1565     }
1566     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1567     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1568     if (!track->PropagateTo(r[0])) {
1569       isOK=kFALSE;
1570     }
1571     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1572   }
1573   //if (!isOK) { delete track; return 0;}
1574   paramIn = *track;
1575   track->AddCovariance(covar);
1576   //
1577   //
1578   for (Int_t i=imin; i<=imax; i++){
1579     AliTPCclusterMI *c=track->GetClusterPointer(i);
1580     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue;
1581     sector = (c->GetDetector()%18);
1582     if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1583       //continue;
1584     }
1585     ncl++;
1586     //if (RejectCluster(c,track)) continue;
1587     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1588     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1589     if (!track->PropagateTo(r[0])) {
1590       isOK=kFALSE;
1591     }
1592     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1593   }
1594   //if (!isOK) { delete track; return 0;}
1595   paramOut=*track;
1596   //
1597   //
1598   //
1599   if (parin) (*parin)=paramIn;
1600   if (parout) (*parout)=paramOut;
1601   delete track;
1602   return ncl;
1603 }
1604
1605
1606
1607 Bool_t AliTPCseed::RefitTrack(AliTPCseed* /*seed*/, Bool_t /*out*/){
1608   //
1609   //
1610   //
1611   return kFALSE;
1612 }
1613
1614
1615
1616
1617
1618
1619 void  AliTPCseed::GetError(AliTPCclusterMI* cluster, AliExternalTrackParam * param,
1620                                   Double_t& erry, Double_t &errz)
1621 {
1622   //
1623   // Get cluster error at given position
1624   //
1625   AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam();
1626   Double_t tany,tanz;
1627   Double_t snp1=param->GetSnp();
1628   tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1629   //
1630   Double_t tgl1=param->GetTgl();
1631   tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1632   //
1633   Int_t padSize = 0;                          // short pads
1634   if (cluster->GetDetector() >= 36) {
1635     padSize = 1;                              // medium pads
1636     if (cluster->GetRow() > 63) padSize = 2; // long pads
1637   }
1638
1639   erry  = clusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany) );
1640   errz  = clusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) );
1641 }
1642
1643
1644 void  AliTPCseed::GetShape(AliTPCclusterMI* cluster, AliExternalTrackParam * param,
1645                                   Double_t& rmsy, Double_t &rmsz)
1646 {
1647   //
1648   // Get cluster error at given position
1649   //
1650   AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam();
1651   Double_t tany,tanz;
1652   Double_t snp1=param->GetSnp();
1653   tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1654   //
1655   Double_t tgl1=param->GetTgl();
1656   tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1657   //
1658   Int_t padSize = 0;                          // short pads
1659   if (cluster->GetDetector() >= 36) {
1660     padSize = 1;                              // medium pads
1661     if (cluster->GetRow() > 63) padSize = 2; // long pads
1662   }
1663
1664   rmsy  = clusterParam->GetRMSQ( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany), TMath::Abs(cluster->GetMax()) );
1665   rmsz  = clusterParam->GetRMSQ( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) ,TMath::Abs(cluster->GetMax()));
1666 }
1667
1668
1669
1670 Double_t AliTPCseed::GetQCorrGeom(Float_t ty, Float_t tz){
1671   //Geoetrical
1672   //ty    - tangent in local y direction
1673   //tz    - tangent
1674   //
1675   Float_t norm=TMath::Sqrt(1+ty*ty+tz*tz);
1676   return norm;
1677 }
1678
1679 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*/){
1680   //
1681   // Q normalization
1682   //
1683   // return value =  Q Normalization factor
1684   // Normalization - 1 - shape factor part for full drift
1685   //                 1 - electron attachment for 0 drift
1686
1687   // Input parameters:
1688   //
1689   // ipad - 0 short pad
1690   //        1 medium pad
1691   //        2 long pad
1692   //
1693   // type - 0 qmax
1694   //      - 1 qtot
1695   //
1696   //z     - z position (-250,250 cm)
1697   //ty    - tangent in local y direction
1698   //tz    - tangent
1699   //
1700
1701   AliTPCClusterParam * paramCl = AliTPCcalibDB::Instance()->GetClusterParam();
1702   AliTPCParam   * paramTPC = AliTPCcalibDB::Instance()->GetParameters();
1703
1704   if (!paramCl) return 1;
1705   //
1706   Double_t dr =  250.-TMath::Abs(z);
1707   Double_t sy =  paramCl->GetRMS0( 0,ipad, dr, TMath::Abs(ty));
1708   Double_t sy0=  paramCl->GetRMS0(0,ipad, 250, 0);
1709   Double_t sz =  paramCl->GetRMS0( 1,ipad, dr, TMath::Abs(tz));
1710   Double_t sz0=  paramCl->GetRMS0(1,ipad, 250, 0);
1711
1712   Double_t sfactorMax = TMath::Sqrt(sy0*sz0/(sy*sz));
1713
1714
1715   Double_t dt = 1000000*(dr/paramTPC->GetDriftV());  //time in microsecond
1716   Double_t attProb = TMath::Exp(-paramTPC->GetAttCoef()*paramTPC->GetOxyCont()*dt);
1717   //
1718   //
1719   if (type==0) return sfactorMax*attProb;
1720
1721   return attProb;
1722
1723
1724 }
1725
1726
1727 /*
1728 //_______________________________________________________________________
1729 Float_t AliTPCseed::GetTPCClustInfo(Int_t nNeighbours, Int_t type, Int_t row0, Int_t row1, TVectorT<float> *returnVec)
1730 {
1731   //
1732   // TPC cluster information
1733   // type 0: get fraction of found/findable clusters with neighbourhood definition
1734   //      1: found clusters
1735   //      2: findable (number of clusters above and below threshold)
1736   //
1737   // definition of findable clusters:
1738   //            a cluster is defined as findable if there is another cluster
1739   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1740   //           effects with a very simple algorithm.
1741   //
1742
1743   const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam
1744   const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
1745   const Float_t kedgey =3.;
1746
1747   Float_t ncl = 0;
1748   Float_t nclBelowThr = 0; // counts number of clusters below threshold
1749
1750   for (Int_t irow=row0; irow<row1; irow++){
1751     AliTPCclusterMI* cluster = GetClusterPointer(irow);
1752
1753     if (!cluster && irow > 1 && irow < 157) {
1754       Bool_t isClBefore = kFALSE;
1755       Bool_t isClAfter  = kFALSE;
1756       for(Int_t ithres = 1; ithres <= nNeighbours; ithres++) {
1757         AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres);
1758         if (clusterBefore) isClBefore = kTRUE;
1759         AliTPCclusterMI * clusterAfter  = GetClusterPointer(irow + ithres);
1760         if (clusterAfter) isClAfter = kTRUE;
1761       }
1762       if (isClBefore && isClAfter) nclBelowThr++;
1763     }
1764     if (!cluster) continue;
1765     //
1766     //
1767     if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
1768     //
1769     AliTPCTrackerPoint * point = GetTrackPoint(irow);
1770     if (point==0) continue;
1771     Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
1772     if (rsigmay > kClusterShapeCut) continue;
1773     //
1774     if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb
1775     ncl++;
1776   }
1777   if(returnVec->GetNoElements != 3){
1778       returnVec->ResizeTo(3);
1779   }
1780   Float_t nclAll = nclBelowThr+ncl;
1781   returnVec(0) = nclAll>0?ncl/nclAll:0;
1782   returnVec(1) = ncl;
1783   returnVec(2) = nclAll;
1784
1785   if(ncl<10)
1786     return 0;
1787   if(type==0)
1788     if(nclAll>0)
1789       return ncl/nclAll;
1790   if(type==1)
1791     return ncl;
1792   if(type==2)
1793     return nclAll;
1794   return 0;
1795 }
1796 */
1797 //_______________________________________________________________________
1798 Int_t AliTPCseed::GetNumberOfClustersIndices() {
1799   Int_t ncls = 0;
1800   for (int i=0; i < 160; i++) {
1801     if ((fIndex[i] & 0x8000) == 0)
1802       ncls++;
1803   }
1804   return ncls;
1805 }
1806
1807 //_______________________________________________________________________
1808 void AliTPCseed::Clear(Option_t*)
1809 {
1810   // formally seed may allocate memory for clusters (althought this should not happen for
1811   // the seeds in the pool). Hence we need this method for fwd. compatibility
1812   if (fClusterOwner) for (int i=160;i--;) {delete fClusterPointer[i]; fClusterPointer[i] = 0;}
1813 }
1814
1815 TObject* AliTPCseed::Clone(const char* /*newname*/) const
1816 {
1817   // temporary override TObject::Clone to avoid crashes in reco
1818   AliTPCseed* src = (AliTPCseed*)this;
1819   AliTPCseed* dst = new AliTPCseed(*src,fClusterOwner);
1820   return dst;
1821 }