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