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