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