TPC/AliTPCseed.cxx

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