TPC/Rec/AliTPCseed.cxx

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