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