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
  29 ClassImp(AliTPCseed)
  30
  31
  32
  33 AliTPCseed::AliTPCseed():
  34   AliTPCtrack(),
  35   fEsd(0x0),
  36   fClusterOwner(kFALSE),
  37   fRow(0),
  38   fSector(-1),
  39   fRelativeSector(-1),
  40   fCurrentSigmaY2(1e10),
  41   fCurrentSigmaZ2(1e10),
  42   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
  43   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
  44   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
  45   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
  46   //
  47   fErrorY2(1e10),
  48   fErrorZ2(1e10),
  49   fCurrentCluster(0x0),
  50   fCurrentClusterIndex1(-1),
  51   fInDead(kFALSE),
  52   fIsSeeding(kFALSE),
  53   fNoCluster(0),
  54   fSort(0),
  55   fBSigned(kFALSE),
  56   fSeedType(0),
  57   fSeed1(-1),
  58   fSeed2(-1),
  59   fMAngular(0),
  60   fCircular(0),
  61   fClusterMap(159),
  62   fSharedMap(159)
  63 {
  64   //
  65   for (Int_t i=0;i<160;i++) SetClusterIndex2(i,-3);
  66   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
  67   for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
  68   for (Int_t i=0;i<AliPID::kSPECIES;i++)   fTPCr[i]=0.2;
  69   for (Int_t i=0;i<4;i++) {
  70     fDEDX[i] = 0.;
  71     fSDEDX[i] = 1e10;
  72     fNCDEDX[i] = 0;
  73   }
  74   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
  75   //  for (Int_t i=0;i<160;i++) fClusterMap[i]=kFALSE;
  76   //for (Int_t i=0;i<160;i++) fSharedMap[i]=kFALSE;
  77   fClusterMap.ResetAllBits(kFALSE);
  78   fSharedMap.ResetAllBits(kFALSE);
  79
  80 }
  81
  82 AliTPCseed::AliTPCseed(const AliTPCseed &s, Bool_t clusterOwner):
  83   AliTPCtrack(s),
  84   fEsd(0x0),
  85   fClusterOwner(clusterOwner),
  86   fRow(0),
  87   fSector(-1),
  88   fRelativeSector(-1),
  89   fCurrentSigmaY2(-1),
  90   fCurrentSigmaZ2(-1),
  91   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
  92   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
  93   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
  94   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
  95   fErrorY2(1e10),
  96   fErrorZ2(1e10),
  97   fCurrentCluster(0x0),
  98   fCurrentClusterIndex1(-1),
  99   fInDead(kFALSE),
 100   fIsSeeding(kFALSE),
 101   fNoCluster(0),
 102   fSort(0),
 103   fBSigned(kFALSE),
 104   fSeedType(0),
 105   fSeed1(-1),
 106   fSeed2(-1),
 107   fMAngular(0),
 108   fCircular(0),
 109   fClusterMap(s.fClusterMap),
 110   fSharedMap(s.fSharedMap)
 111 {
 112   //---------------------
 113   // dummy copy constructor
 114   //-------------------------
 115   for (Int_t i=0;i<160;i++) {
 116     fClusterPointer[i]=0;
 117     if (fClusterOwner){
 118       if (s.fClusterPointer[i])
 119         fClusterPointer[i] = new AliTPCclusterMI(*(s.fClusterPointer[i]));
 120     }else{
 121       fClusterPointer[i] = s.fClusterPointer[i];
 122     }
 123     fTrackPoints[i] = s.fTrackPoints[i];
 124   }
 125   for (Int_t i=0;i<160;i++) fIndex[i] = s.fIndex[i];
 126   for (Int_t i=0;i<AliPID::kSPECIES;i++)   fTPCr[i]=s.fTPCr[i];
 127   for (Int_t i=0;i<4;i++) {
 128     fDEDX[i] = s.fDEDX[i];
 129     fSDEDX[i] = s.fSDEDX[i];
 130     fNCDEDX[i] = s.fNCDEDX[i];
 131   }
 132   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = s.fOverlapLabels[i];
 133
 134 }
 135
 136
 137 AliTPCseed::AliTPCseed(const AliTPCtrack &t):
 138   AliTPCtrack(t),
 139   fEsd(0x0),
 140   fClusterOwner(kFALSE),
 141   fRow(0),
 142   fSector(-1),
 143   fRelativeSector(-1),
 144   fCurrentSigmaY2(-1),
 145   fCurrentSigmaZ2(-1),
 146   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
 147   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
 148   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
 149   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
 150   fErrorY2(1e10),
 151   fErrorZ2(1e10),
 152   fCurrentCluster(0x0),
 153   fCurrentClusterIndex1(-1),
 154   fInDead(kFALSE),
 155   fIsSeeding(kFALSE),
 156   fNoCluster(0),
 157   fSort(0),
 158   fBSigned(kFALSE),
 159   fSeedType(0),
 160   fSeed1(-1),
 161   fSeed2(-1),
 162   fMAngular(0),
 163   fCircular(0),
 164   fClusterMap(159),
 165   fSharedMap(159)
 166 {
 167   //
 168   // Constructor from AliTPCtrack
 169   //
 170   fFirstPoint =0;
 171   for (Int_t i=0;i<5;i++)   fTPCr[i]=0.2;
 172   for (Int_t i=0;i<160;i++) {
 173     fClusterPointer[i] = 0;
 174     Int_t index = t.GetClusterIndex(i);
 175     if (index>=-1){
 176       SetClusterIndex2(i,index);
 177     }
 178     else{
 179       SetClusterIndex2(i,-3);
 180     }
 181   }
 182   for (Int_t i=0;i<4;i++) {
 183     fDEDX[i] = 0.;
 184     fSDEDX[i] = 1e10;
 185     fNCDEDX[i] = 0;
 186   }
 187   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
 188
 189   //for (Int_t i=0;i<160;i++) fClusterMap[i]=kFALSE;
 190   //for (Int_t i=0;i<160;i++) fSharedMap[i]=kFALSE;
 191   fClusterMap.ResetAllBits(kFALSE);
 192   fSharedMap.ResetAllBits(kFALSE);
 193
 194 }
 195
 196 AliTPCseed::AliTPCseed(Double_t xr, Double_t alpha, const Double_t xx[5],
 197                        const Double_t cc[15], Int_t index):
 198   AliTPCtrack(xr, alpha, xx, cc, index),
 199   fEsd(0x0),
 200   fClusterOwner(kFALSE),
 201   fRow(0),
 202   fSector(-1),
 203   fRelativeSector(-1),
 204   fCurrentSigmaY2(-1),
 205   fCurrentSigmaZ2(-1),
 206   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
 207   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
 208   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
 209   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
 210   fErrorY2(1e10),
 211   fErrorZ2(1e10),
 212   fCurrentCluster(0x0),
 213   fCurrentClusterIndex1(-1),
 214   fInDead(kFALSE),
 215   fIsSeeding(kFALSE),
 216   fNoCluster(0),
 217   fSort(0),
 218   fBSigned(kFALSE),
 219   fSeedType(0),
 220   fSeed1(-1),
 221   fSeed2(-1),
 222   fMAngular(0),
 223   fCircular(0),
 224   fClusterMap(159),
 225   fSharedMap(159)
 226 {
 227   //
 228   // Constructor
 229   //
 230   fFirstPoint =0;
 231   for (Int_t i=0;i<160;i++) SetClusterIndex2(i,-3);
 232   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
 233   for (Int_t i=0;i<5;i++)   fTPCr[i]=0.2;
 234   for (Int_t i=0;i<4;i++) {
 235     fDEDX[i] = 0.;
 236     fSDEDX[i] = 1e10;
 237     fNCDEDX[i] = 0;
 238   }
 239   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
 240 }
 241
 242 AliTPCseed::~AliTPCseed(){
 243   //
 244   // destructor
 245   fNoCluster =0;
 246   if (fClusterOwner){
 247     for (Int_t icluster=0; icluster<160; icluster++){
 248       delete fClusterPointer[icluster];
 249     }
 250   }
 251
 252 }
 253 //_________________________________________________
 254 AliTPCseed & AliTPCseed::operator=(const AliTPCseed &param)
 255 {
 256   //
 257   // assignment operator
 258   //
 259   if(this!=&param){
 260     AliTPCtrack::operator=(param);
 261     fEsd =param.fEsd;
 262     for(Int_t i = 0;i<160;++i)fClusterPointer[i] = param.fClusterPointer[i]; // this is not allocated by AliTPCSeed
 263     fClusterOwner = param.fClusterOwner;
 264     // leave out fPoint, they are also not copied in the copy ctor...
 265     // but deleted in the dtor... strange...
 266     fRow            = param.fRow;
 267     fSector         = param.fSector;
 268     fRelativeSector = param.fRelativeSector;
 269     fCurrentSigmaY2 = param.fCurrentSigmaY2;
 270     fCurrentSigmaZ2 = param.fCurrentSigmaZ2;
 271     fErrorY2        = param.fErrorY2;
 272     fErrorZ2        = param.fErrorZ2;
 273     fCurrentCluster = param.fCurrentCluster; // this is not allocated by AliTPCSeed
 274     fCurrentClusterIndex1 = param.fCurrentClusterIndex1;
 275     fInDead         = param.fInDead;
 276     fIsSeeding      = param.fIsSeeding;
 277     fNoCluster      = param.fNoCluster;
 278     fSort           = param.fSort;
 279     fBSigned        = param.fBSigned;
 280     for(Int_t i = 0;i<4;++i){
 281       fDEDX[i]   = param.fDEDX[i];
 282       fSDEDX[i]  = param.fSDEDX[i];
 283       fNCDEDX[i] = param.fNCDEDX[i];
 284     }
 285     for(Int_t i = 0;i<AliPID::kSPECIES;++i)fTPCr[i] = param.fTPCr[i];
 286
 287     fSeedType = param.fSeedType;
 288     fSeed1    = param.fSeed1;
 289     fSeed2    = param.fSeed2;
 290     for(Int_t i = 0;i<12;++i)fOverlapLabels[i] = param.fOverlapLabels[i];
 291     fMAngular = param.fMAngular;
 292     fCircular = param.fCircular;
 293     for(int i = 0;i<160;++i)fTrackPoints[i] =  param.fTrackPoints[i];
 294     fClusterMap = param.fClusterMap;
 295     fSharedMap = param.fSharedMap;
 296   }
 297   return (*this);
 298 }
 299 //____________________________________________________
 300 AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
 301 {
 302   //
 303   //
 304   return &fTrackPoints[i];
 305 }
 306
 307
 308
 309 Double_t AliTPCseed::GetDensityFirst(Int_t n)
 310 {
 311   //
 312   //
 313   // return cluster for n rows bellow first point
 314   Int_t nfoundable = 1;
 315   Int_t nfound      = 1;
 316   for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
 317     Int_t index = GetClusterIndex2(i);
 318     if (index!=-1) nfoundable++;
 319     if (index>0) nfound++;
 320   }
 321   if (nfoundable<n) return 0;
 322   return Double_t(nfound)/Double_t(nfoundable);
 323
 324 }
 325
 326
 327 void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
 328 {
 329   // get cluster stat.  on given region
 330   //
 331   found       = 0;
 332   foundable   = 0;
 333   shared      =0;
 334   for (Int_t i=first;i<last; i++){
 335     Int_t index = GetClusterIndex2(i);
 336     if (index!=-1) foundable++;
 337     if (index&0x8000) continue;
 338     if (fClusterPointer[i]) {
 339       found++;
 340     }
 341     else
 342       continue;
 343
 344     if (fClusterPointer[i]->IsUsed(10)) {
 345       shared++;
 346       continue;
 347     }
 348     if (!plus2) continue; //take also neighborhoud
 349     //
 350     if ( (i>0) && fClusterPointer[i-1]){
 351       if (fClusterPointer[i-1]->IsUsed(10)) {
 352         shared++;
 353         continue;
 354       }
 355     }
 356     if ( fClusterPointer[i+1]){
 357       if (fClusterPointer[i+1]->IsUsed(10)) {
 358         shared++;
 359         continue;
 360       }
 361     }
 362
 363   }
 364   //if (shared>found){
 365     //Error("AliTPCseed::GetClusterStatistic","problem\n");
 366   //}
 367 }
 368
 369
 370
 371
 372
 373 void AliTPCseed::Reset(Bool_t all)
 374 {
 375   //
 376   //
 377   SetNumberOfClusters(0);
 378   fNFoundable = 0;
 379   SetChi2(0);
 380   ResetCovariance(10.);
 381   /*
 382   if (fTrackPoints){
 383     for (Int_t i=0;i<8;i++){
 384       delete [] fTrackPoints[i];
 385     }
 386     delete fTrackPoints;
 387     fTrackPoints =0;
 388   }
 389   */
 390
 391   if (all){
 392     for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
 393     for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
 394   }
 395
 396 }
 397
 398
 399 void AliTPCseed::Modify(Double_t factor)
 400 {
 401
 402   //------------------------------------------------------------------
 403   //This function makes a track forget its history :)
 404   //------------------------------------------------------------------
 405   if (factor<=0) {
 406     ResetCovariance(10.);
 407     return;
 408   }
 409   ResetCovariance(factor);
 410
 411   SetNumberOfClusters(0);
 412   fNFoundable =0;
 413   SetChi2(0);
 414   fRemoval = 0;
 415   fCurrentSigmaY2 = 0.000005;
 416   fCurrentSigmaZ2 = 0.000005;
 417   fNoCluster     = 0;
 418   //fFirstPoint = 160;
 419   //fLastPoint  = 0;
 420 }
 421
 422
 423
 424
 425 Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
 426 {
 427   //-----------------------------------------------------------------
 428   // This function find proloncation of a track to a reference plane x=xk.
 429   // doesn't change internal state of the track
 430   //-----------------------------------------------------------------
 431
 432   Double_t x1=GetX(), x2=x1+(xk-x1), dx=x2-x1;
 433
 434   if (TMath::Abs(GetSnp()+GetC()*dx) >= AliTPCReconstructor::GetMaxSnpTrack()) {
 435     return 0;
 436   }
 437
 438   //  Double_t y1=fP0, z1=fP1;
 439   Double_t c1=GetSnp(), r1=sqrt(1.- c1*c1);
 440   Double_t c2=c1 + GetC()*dx, r2=sqrt(1.- c2*c2);
 441
 442   y = GetY();
 443   z = GetZ();
 444   //y += dx*(c1+c2)/(r1+r2);
 445   //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
 446
 447   Double_t dy = dx*(c1+c2)/(r1+r2);
 448   Double_t dz = 0;
 449   //
 450   Double_t delta = GetC()*dx*(c1+c2)/(c1*r2 + c2*r1);
 451   /*
 452   if (TMath::Abs(delta)>0.0001){
 453     dz = fP3*TMath::ASin(delta)/fP4;
 454   }else{
 455     dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
 456   }
 457   */
 458   //  dz =  fP3*AliTPCFastMath::FastAsin(delta)/fP4;
 459   dz =  GetTgl()*TMath::ASin(delta)/GetC();
 460   //
 461   y+=dy;
 462   z+=dz;
 463
 464
 465   return 1;
 466 }
 467
 468
 469 //_____________________________________________________________________________
 470 Double_t AliTPCseed::GetPredictedChi2(const AliCluster *c) const
 471 {
 472   //-----------------------------------------------------------------
 473   // This function calculates a predicted chi2 increment.
 474   //-----------------------------------------------------------------
 475   Double_t p[2]={c->GetY(), c->GetZ()};
 476   Double_t cov[3]={fErrorY2, 0., fErrorZ2};
 477   return AliExternalTrackParam::GetPredictedChi2(p,cov);
 478 }
 479
 480 //_________________________________________________________________________________________
 481
 482
 483 Int_t AliTPCseed::Compare(const TObject *o) const {
 484   //-----------------------------------------------------------------
 485   // This function compares tracks according to the sector - for given sector according z
 486   //-----------------------------------------------------------------
 487   AliTPCseed *t=(AliTPCseed*)o;
 488
 489   if (fSort == 0){
 490     if (t->fRelativeSector>fRelativeSector) return -1;
 491     if (t->fRelativeSector<fRelativeSector) return 1;
 492     Double_t z2 = t->GetZ();
 493     Double_t z1 = GetZ();
 494     if (z2>z1) return 1;
 495     if (z2<z1) return -1;
 496     return 0;
 497   }
 498   else {
 499     Float_t f2 =1;
 500     f2 = 1-20*TMath::Sqrt(t->GetSigma1Pt2())/(t->OneOverPt()+0.0066);
 501     if (t->fBConstrain) f2=1.2;
 502
 503     Float_t f1 =1;
 504     f1 = 1-20*TMath::Sqrt(GetSigma1Pt2())/(OneOverPt()+0.0066);
 505
 506     if (fBConstrain)   f1=1.2;
 507
 508     if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
 509     else return +1;
 510   }
 511 }
 512
 513
 514
 515
 516 //_____________________________________________________________________________
 517 Bool_t AliTPCseed::Update(const AliCluster *c, Double_t chisq, Int_t index)
 518 {
 519   //-----------------------------------------------------------------
 520   // This function associates a cluster with this track.
 521   //-----------------------------------------------------------------
 522   Int_t n=GetNumberOfClusters();
 523   Int_t idx=GetClusterIndex(n);    // save the current cluster index
 524
 525   AliCluster cl(*c);  cl.SetSigmaY2(fErrorY2); cl.SetSigmaZ2(fErrorZ2);
 526   if (!AliTPCtrack::Update(&cl,chisq,index)) return kFALSE;
 527
 528   if (fCMeanSigmaY2p30<0){
 529     fCMeanSigmaY2p30= c->GetSigmaY2();   //! current mean sigma Y2 - mean30%
 530     fCMeanSigmaZ2p30= c->GetSigmaZ2();   //! current mean sigma Z2 - mean30%
 531     fCMeanSigmaY2p30R = 1;   //! current mean sigma Y2 - mean5%
 532     fCMeanSigmaZ2p30R = 1;   //! current mean sigma Z2 - mean5%
 533   }
 534   //
 535   fCMeanSigmaY2p30= 0.70*fCMeanSigmaY2p30 +0.30*c->GetSigmaY2();
 536   fCMeanSigmaZ2p30= 0.70*fCMeanSigmaZ2p30 +0.30*c->GetSigmaZ2();
 537   if (fCurrentSigmaY2>0){
 538     fCMeanSigmaY2p30R = 0.7*fCMeanSigmaY2p30R  +0.3*c->GetSigmaY2()/fCurrentSigmaY2;
 539     fCMeanSigmaZ2p30R = 0.7*fCMeanSigmaZ2p30R  +0.3*c->GetSigmaZ2()/fCurrentSigmaZ2;
 540   }
 541
 542
 543   SetClusterIndex(n,idx);          // restore the current cluster index
 544   return kTRUE;
 545 }
 546
 547
 548
 549 //_____________________________________________________________________________
 550 Float_t AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
 551   //-----------------------------------------------------------------
 552   // This funtion calculates dE/dX within the "low" and "up" cuts.
 553   //-----------------------------------------------------------------
 554
 555   Float_t amp[200];
 556   Float_t angular[200];
 557   Float_t weight[200];
 558   Int_t index[200];
 559   //Int_t nc = 0;
 560   Float_t meanlog = 100.;
 561
 562   Float_t mean[4]  = {0,0,0,0};
 563   Float_t sigma[4] = {1000,1000,1000,1000};
 564   Int_t nc[4]      = {0,0,0,0};
 565   Float_t norm[4]    = {1000,1000,1000,1000};
 566   //
 567   //
 568   fNShared =0;
 569
 570   for (Int_t of =0; of<4; of++){
 571     for (Int_t i=of+i1;i<i2;i+=4)
 572       {
 573         Int_t index = fIndex[i];
 574         if (index<0||index&0x8000) continue;
 575
 576         //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
 577         AliTPCTrackerPoint * point = GetTrackPoint(i);
 578         //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
 579         //AliTPCTrackerPoint * pointp = 0;
 580         //if (i<159) pointp = GetTrackPoint(i+1);
 581
 582         if (point==0) continue;
 583         AliTPCclusterMI * cl = fClusterPointer[i];
 584         if (cl==0) continue;
 585         if (onlyused && (!cl->IsUsed(10))) continue;
 586         if (cl->IsUsed(11)) {
 587           fNShared++;
 588           continue;
 589         }
 590         Int_t   type   = cl->GetType();
 591         //if (point->fIsShared){
 592         //  fNShared++;
 593         //  continue;
 594         //}
 595         //if (pointm)
 596         //  if (pointm->fIsShared) continue;
 597         //if (pointp)
 598         //  if (pointp->fIsShared) continue;
 599
 600         if (type<0) continue;
 601         //if (type>10) continue;
 602         //if (point->GetErrY()==0) continue;
 603         //if (point->GetErrZ()==0) continue;
 604
 605         //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
 606         //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
 607         //if ((ddy*ddy+ddz*ddz)>10) continue;
 608
 609
 610         //      if (point->GetCPoint().GetMax()<5) continue;
 611         if (cl->GetMax()<5) continue;
 612         Float_t angley = point->GetAngleY();
 613         Float_t anglez = point->GetAngleZ();
 614
 615         Float_t rsigmay2 =  point->GetSigmaY();
 616         Float_t rsigmaz2 =  point->GetSigmaZ();
 617         /*
 618         Float_t ns = 1.;
 619         if (pointm){
 620           rsigmay +=  pointm->GetTPoint().GetSigmaY();
 621           rsigmaz +=  pointm->GetTPoint().GetSigmaZ();
 622           ns+=1.;
 623         }
 624         if (pointp){
 625           rsigmay +=  pointp->GetTPoint().GetSigmaY();
 626           rsigmaz +=  pointp->GetTPoint().GetSigmaZ();
 627           ns+=1.;
 628         }
 629         rsigmay/=ns;
 630         rsigmaz/=ns;
 631         */
 632
 633         Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
 634
 635         Float_t ampc   = 0;     // normalization to the number of electrons
 636         if (i>64){
 637           //      ampc = 1.*point->GetCPoint().GetMax();
 638           ampc = 1.*cl->GetMax();
 639           //ampc = 1.*point->GetCPoint().GetQ();
 640           //      AliTPCClusterPoint & p = point->GetCPoint();
 641           //      Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
 642           // Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
 643           //Float_t dz =
 644           //  TMath::Abs( Int_t(iz) - iz + 0.5);
 645           //ampc *= 1.15*(1-0.3*dy);
 646           //ampc *= 1.15*(1-0.3*dz);
 647           //      Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
 648           //ampc               *=zfactor;
 649         }
 650         else{
 651           //ampc = 1.0*point->GetCPoint().GetMax();
 652           ampc = 1.0*cl->GetMax();
 653           //ampc = 1.0*point->GetCPoint().GetQ();
 654           //AliTPCClusterPoint & p = point->GetCPoint();
 655           // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
 656           //Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
 657           //Float_t dz =
 658           //  TMath::Abs( Int_t(iz) - iz + 0.5);
 659
 660           //ampc *= 1.15*(1-0.3*dy);
 661           //ampc *= 1.15*(1-0.3*dz);
 662           //    Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
 663           //ampc               *=zfactor;
 664
 665         }
 666         ampc *= 2.0;     // put mean value to channel 50
 667         //ampc *= 0.58;     // put mean value to channel 50
 668         Float_t w      =  1.;
 669         //      if (type>0)  w =  1./(type/2.-0.5);
 670         //      Float_t z = TMath::Abs(cl->GetZ());
 671         if (i<64) {
 672           ampc /= 0.6;
 673           //ampc /= (1+0.0008*z);
 674         } else
 675           if (i>128){
 676             ampc /=1.5;
 677             //ampc /= (1+0.0008*z);
 678           }else{
 679             //ampc /= (1+0.0008*z);
 680           }
 681
 682         if (type<0) {  //amp at the border - lower weight
 683           // w*= 2.;
 684
 685           continue;
 686         }
 687         if (rsigma>1.5) ampc/=1.3;  // if big backround
 688         amp[nc[of]]        = ampc;
 689         angular[nc[of]]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
 690         weight[nc[of]]     = w;
 691         nc[of]++;
 692       }
 693
 694     TMath::Sort(nc[of],amp,index,kFALSE);
 695     Float_t sumamp=0;
 696     Float_t sumamp2=0;
 697     Float_t sumw=0;
 698     //meanlog = amp[index[Int_t(nc[of]*0.33)]];
 699     meanlog = 50;
 700     for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
 701       Float_t ampl      = amp[index[i]]/angular[index[i]];
 702       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
 703       //
 704       sumw    += weight[index[i]];
 705       sumamp  += weight[index[i]]*ampl;
 706       sumamp2 += weight[index[i]]*ampl*ampl;
 707       norm[of]    += angular[index[i]]*weight[index[i]];
 708     }
 709     if (sumw<1){
 710       SetdEdx(0);
 711     }
 712     else {
 713       norm[of] /= sumw;
 714       mean[of]  = sumamp/sumw;
 715       sigma[of] = sumamp2/sumw-mean[of]*mean[of];
 716       if (sigma[of]>0.1)
 717         sigma[of] = TMath::Sqrt(sigma[of]);
 718       else
 719         sigma[of] = 1000;
 720
 721     mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
 722     //mean  *=(1-0.02*(sigma/(mean*0.17)-1.));
 723     //mean *=(1-0.1*(norm-1.));
 724     }
 725   }
 726
 727   Float_t dedx =0;
 728   fSdEdx =0;
 729   fMAngular =0;
 730   //  mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
 731   //  mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
 732
 733
 734   //  dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
 735   //  (  TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
 736
 737   Int_t norm2 = 0;
 738   Int_t norm3 = 0;
 739   for (Int_t i =0;i<4;i++){
 740     if (nc[i]>2&&nc[i]<1000){
 741       dedx      += mean[i] *nc[i];
 742       fSdEdx    += sigma[i]*(nc[i]-2);
 743       fMAngular += norm[i] *nc[i];
 744       norm2     += nc[i];
 745       norm3     += nc[i]-2;
 746     }
 747     fDEDX[i]  = mean[i];
 748     fSDEDX[i] = sigma[i];
 749     fNCDEDX[i]= nc[i];
 750   }
 751
 752   if (norm3>0){
 753     dedx   /=norm2;
 754     fSdEdx /=norm3;
 755     fMAngular/=norm2;
 756   }
 757   else{
 758     SetdEdx(0);
 759     return 0;
 760   }
 761   //  Float_t dedx1 =dedx;
 762   /*
 763   dedx =0;
 764   for (Int_t i =0;i<4;i++){
 765     if (nc[i]>2&&nc[i]<1000){
 766       mean[i]   = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
 767       dedx      += mean[i] *nc[i];
 768     }
 769     fDEDX[i]  = mean[i];
 770   }
 771   dedx /= norm2;
 772   */
 773
 774
 775   SetdEdx(dedx);
 776   return dedx;
 777 }
 778 Double_t AliTPCseed::Bethe(Double_t bg){
 779   //
 780   // This is the Bethe-Bloch function normalised to 1 at the minimum
 781   //
 782   Double_t bg2=bg*bg;
 783   Double_t bethe;
 784   if (bg<3.5e1)
 785     bethe=(1.+ bg2)/bg2*(log(5940*bg2) - bg2/(1.+ bg2));
 786   else // Density effect ( approximately :)
 787     bethe=1.15*(1.+ bg2)/bg2*(log(3.5*5940*bg) - bg2/(1.+ bg2));
 788   return bethe/11.091;
 789 }
 790
 791 void AliTPCseed::CookPID()
 792 {
 793   //
 794   // cook PID information according dEdx
 795   //
 796   Double_t fRange = 10.;
 797   Double_t fRes   = 0.1;
 798   Double_t fMIP   = 47.;
 799   //
 800   Int_t ns=AliPID::kSPECIES;
 801   Double_t sumr =0;
 802   for (Int_t j=0; j<ns; j++) {
 803     Double_t mass=AliPID::ParticleMass(j);
 804     Double_t mom=GetP();
 805     Double_t dedx=fdEdx/fMIP;
 806     Double_t bethe=Bethe(mom/mass);
 807     Double_t sigma=fRes*bethe;
 808     if (sigma>0.001){
 809       if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 810         fTPCr[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 811         sumr+=fTPCr[j];
 812         continue;
 813       }
 814       fTPCr[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 815       sumr+=fTPCr[j];
 816     }
 817     else{
 818       fTPCr[j]=1.;
 819       sumr+=fTPCr[j];
 820     }
 821   }
 822   for (Int_t j=0; j<ns; j++) {
 823     fTPCr[j]/=sumr;           //normalize
 824   }
 825 }
 826
 827 /*
 828 void AliTPCseed::CookdEdx2(Double_t low, Double_t up) {
 829   //-----------------------------------------------------------------
 830   // This funtion calculates dE/dX within the "low" and "up" cuts.
 831   //-----------------------------------------------------------------
 832
 833   Float_t amp[200];
 834   Float_t angular[200];
 835   Float_t weight[200];
 836   Int_t index[200];
 837   Bool_t inlimit[200];
 838   for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
 839   for (Int_t i=0;i<200;i++) amp[i]=10000;
 840   for (Int_t i=0;i<200;i++) angular[i]= 1;;
 841
 842
 843   //
 844   Float_t meanlog = 100.;
 845   Int_t indexde[4]={0,64,128,160};
 846
 847   Float_t amean     =0;
 848   Float_t asigma    =0;
 849   Float_t anc       =0;
 850   Float_t anorm     =0;
 851
 852   Float_t mean[4]  = {0,0,0,0};
 853   Float_t sigma[4] = {1000,1000,1000,1000};
 854   Int_t nc[4]      = {0,0,0,0};
 855   Float_t norm[4]    = {1000,1000,1000,1000};
 856   //
 857   //
 858   fNShared =0;
 859
 860   //  for (Int_t of =0; of<3; of++){
 861   //  for (Int_t i=indexde[of];i<indexde[of+1];i++)
 862   for (Int_t i =0; i<160;i++)
 863     {
 864         AliTPCTrackPoint * point = GetTrackPoint(i);
 865         if (point==0) continue;
 866         if (point->fIsShared){
 867           fNShared++;
 868           continue;
 869         }
 870         Int_t   type   = point->GetCPoint().GetType();
 871         if (type<0) continue;
 872         if (point->GetCPoint().GetMax()<5) continue;
 873         Float_t angley = point->GetTPoint().GetAngleY();
 874         Float_t anglez = point->GetTPoint().GetAngleZ();
 875         Float_t rsigmay =  point->GetCPoint().GetSigmaY();
 876         Float_t rsigmaz =  point->GetCPoint().GetSigmaZ();
 877         Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
 878
 879         Float_t ampc   = 0;     // normalization to the number of electrons
 880         if (i>64){
 881           ampc =  point->GetCPoint().GetMax();
 882         }
 883         else{
 884           ampc = point->GetCPoint().GetMax();
 885         }
 886         ampc *= 2.0;     // put mean value to channel 50
 887         //      ampc *= 0.565;     // put mean value to channel 50
 888
 889         Float_t w      =  1.;
 890         Float_t z = TMath::Abs(point->GetCPoint().GetZ());
 891         if (i<64) {
 892           ampc /= 0.63;
 893         } else
 894           if (i>128){
 895             ampc /=1.51;
 896           }
 897         if (type<0) {  //amp at the border - lower weight
 898           continue;
 899         }
 900         if (rsigma>1.5) ampc/=1.3;  // if big backround
 901         angular[i]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
 902         amp[i]        = ampc/angular[i];
 903         weight[i]     = w;
 904         anc++;
 905     }
 906
 907   TMath::Sort(159,amp,index,kFALSE);
 908   for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){
 909     inlimit[index[i]] = kTRUE;  // take all clusters
 910   }
 911
 912   //  meanlog = amp[index[Int_t(anc*0.3)]];
 913   meanlog =10000.;
 914   for (Int_t of =0; of<3; of++){
 915     Float_t sumamp=0;
 916     Float_t sumamp2=0;
 917     Float_t sumw=0;
 918    for (Int_t i=indexde[of];i<indexde[of+1];i++)
 919       {
 920         if (inlimit[i]==kFALSE) continue;
 921         Float_t ampl      = amp[i];
 922         ///angular[i];
 923         ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
 924         //
 925         sumw    += weight[i];
 926         sumamp  += weight[i]*ampl;
 927         sumamp2 += weight[i]*ampl*ampl;
 928         norm[of]    += angular[i]*weight[i];
 929         nc[of]++;
 930       }
 931    if (sumw<1){
 932      SetdEdx(0);
 933    }
 934    else {
 935      norm[of] /= sumw;
 936      mean[of]  = sumamp/sumw;
 937      sigma[of] = sumamp2/sumw-mean[of]*mean[of];
 938      if (sigma[of]>0.1)
 939        sigma[of] = TMath::Sqrt(sigma[of]);
 940      else
 941        sigma[of] = 1000;
 942      mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
 943    }
 944   }
 945
 946   Float_t dedx =0;
 947   fSdEdx =0;
 948   fMAngular =0;
 949   //
 950   Int_t norm2 = 0;
 951   Int_t norm3 = 0;
 952   Float_t www[3] = {12.,14.,17.};
 953   //Float_t www[3] = {1.,1.,1.};
 954
 955   for (Int_t i =0;i<3;i++){
 956     if (nc[i]>2&&nc[i]<1000){
 957       dedx      += mean[i] *nc[i]*www[i]/sigma[i];
 958       fSdEdx    += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
 959       fMAngular += norm[i] *nc[i];
 960       norm2     += nc[i]*www[i]/sigma[i];
 961       norm3     += (nc[i]-2)*www[i]/sigma[i];
 962     }
 963     fDEDX[i]  = mean[i];
 964     fSDEDX[i] = sigma[i];
 965     fNCDEDX[i]= nc[i];
 966   }
 967
 968   if (norm3>0){
 969     dedx   /=norm2;
 970     fSdEdx /=norm3;
 971     fMAngular/=norm2;
 972   }
 973   else{
 974     SetdEdx(0);
 975     return;
 976   }
 977   //  Float_t dedx1 =dedx;
 978
 979   dedx =0;
 980   Float_t norm4 = 0;
 981   for (Int_t i =0;i<3;i++){
 982     if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
 983       //mean[i]   = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
 984       dedx      += mean[i] *(nc[i])/(sigma[i]);
 985       norm4     += (nc[i])/(sigma[i]);
 986     }
 987     fDEDX[i]  = mean[i];
 988   }
 989   if (norm4>0) dedx /= norm4;
 990
 991
 992
 993   SetdEdx(dedx);
 994
 995   //mi deDX
 996
 997 }
 998 */
 999 Double_t AliTPCseed::GetYat(Double_t xk) const {
1000 //-----------------------------------------------------------------
1001 // This function calculates the Y-coordinate of a track at the plane x=xk.
1002 //-----------------------------------------------------------------
1003   if (TMath::Abs(GetSnp())>AliTPCReconstructor::GetMaxSnpTrack()) return 0.; //patch 01 jan 06
1004     Double_t c1=GetSnp(), r1=TMath::Sqrt(1.- c1*c1);
1005     Double_t c2=c1+GetC()*(xk-GetX());
1006     if (TMath::Abs(c2)>AliTPCReconstructor::GetMaxSnpTrack()) return 0;
1007     Double_t r2=TMath::Sqrt(1.- c2*c2);
1008     return GetY() + (xk-GetX())*(c1+c2)/(r1+r2);
1009 }
1010
1011 void AliTPCseed::SetClusterMapBit(int ibit, Bool_t state)
1012 {
1013   fClusterMap[ibit] = state;
1014 }
1015 Bool_t AliTPCseed::GetClusterMapBit(int ibit)
1016 {
1017   return fClusterMap[ibit];
1018 }
1019 void AliTPCseed::SetSharedMapBit(int ibit, Bool_t state)
1020 {
1021   fSharedMap[ibit] = state;
1022 }
1023 Bool_t AliTPCseed::GetSharedMapBit(int ibit)
1024 {
1025   return fSharedMap[ibit];
1026 }
1027
1028
1029
1030
1031
1032 Float_t  AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2){
1033
1034   //
1035   // calculates dedx using the cluster
1036   // low    -  up specify trunc mean range  - default form 0-0.7
1037   // type   -  1 - max charge  or 0- total charge in cluster
1038   //           //2- max no corr 3- total+ correction
1039   // i1-i2  -  the pad-row range used for calculation
1040   //
1041   // normalization parametrization taken from AliTPCClusterParam
1042   //
1043   AliTPCClusterParam * parcl = AliTPCClusterParam::Instance();
1044   if (!parcl) return 0;
1045   Float_t amp[160];
1046   Int_t   indexes[160];
1047   Int_t   ncl=0;
1048   //
1049   //
1050   const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
1051   const Float_t kedgey =4.;
1052   //
1053   for (Int_t irow=i1; irow<i2; irow++){
1054     AliTPCclusterMI* cluster = GetClusterPointer(irow);
1055     if (!cluster) continue;
1056     if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
1057     Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ();
1058     //do normalization
1059     Float_t corr=1;
1060     Int_t  ipad= 0;
1061     if (irow>62) ipad=1;
1062     if (irow>127) ipad=2;
1063     if (type<=1){
1064       //
1065       AliTPCTrackerPoint * point = GetTrackPoint(irow);
1066       Float_t              ty = TMath::Abs(point->GetAngleY());
1067       Float_t              tz = TMath::Abs(point->GetAngleZ());
1068
1069       Float_t dr    = (250.-TMath::Abs(cluster->GetZ()))/250.;
1070       corr  = parcl->Qnorm(ipad,type,dr,ty,tz);
1071     }
1072     amp[ncl]=charge/corr;
1073
1074     amp[ncl] *= 2.0;     // put mean value to channel 50
1075     if (ipad==0) {
1076       amp[ncl] /= 0.65; // this we will take form OCDB
1077     } else
1078       if (ipad==2){
1079         amp[ncl] /=1.57;
1080       }else{
1081       }
1082     ncl++;
1083   }
1084
1085   if (type>3) return ncl;
1086   TMath::Sort(ncl,amp, indexes, kFALSE);
1087
1088   if (ncl<10) return 0;
1089
1090   Float_t suma=0;
1091   Float_t sumn=0;
1092   Int_t icl0=TMath::Nint(ncl*low);
1093   Int_t icl1=TMath::Nint(ncl*up);
1094   for (Int_t icl=icl0; icl<icl1;icl++){
1095     suma+=amp[indexes[icl]];
1096     sumn++;
1097   }
1098   return suma/sumn;
1099
1100 }
1101
1102 Double_t AliTPCseed::BetheMass(Double_t mass){
1103   //
1104   // return bethe-bloch
1105   //
1106   Float_t bg= P()/mass;
1107   const Double_t kp1=0.76176e-1;
1108   const Double_t kp2=10.632;
1109   const Double_t kp3=0.13279e-4;
1110   const Double_t kp4=1.8631;
1111   const Double_t kp5=1.9479;
1112
1113   Double_t dbg = (Double_t) bg;
1114
1115   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
1116
1117   Double_t aa = TMath::Power(beta,kp4);
1118   Double_t bb = TMath::Power(1./dbg,kp5);
1119
1120   bb=TMath::Log(kp3+bb);
1121
1122   return ((Float_t)((kp2-aa-bb)*kp1/aa));
1123 }
1124
1125
1126 Float_t  AliTPCseed::CookShape(Int_t type){
1127   //
1128   //
1129   //
1130  //-----------------------------------------------------------------
1131   // This funtion calculates dE/dX within the "low" and "up" cuts.
1132   //-----------------------------------------------------------------
1133   Float_t means=0;
1134   Float_t meanc=0;
1135   for (Int_t i =0; i<160;i++)    {
1136     AliTPCTrackerPoint * point = GetTrackPoint(i);
1137     if (point==0) continue;
1138
1139     AliTPCclusterMI * cl = fClusterPointer[i];
1140     if (cl==0) continue;
1141
1142     Float_t rsigmay =  TMath::Sqrt(point->GetSigmaY());
1143     Float_t rsigmaz =  TMath::Sqrt(point->GetSigmaZ());
1144     Float_t rsigma =   (rsigmay+rsigmaz)*0.5;
1145     if (type==0) means+=rsigma;
1146     if (type==1) means+=rsigmay;
1147     if (type==2) means+=rsigmaz;
1148     meanc++;
1149   }
1150   Float_t mean = (meanc>0)? means/meanc:0;
1151   return mean;
1152 }