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
  28 ClassImp(AliTPCseed)
  29
  30
  31
  32 AliTPCseed::AliTPCseed():
  33   AliTPCtrack(),
  34   fEsd(0x0),
  35   fPoints(0x0),
  36   fEPoints(0x0),
  37   fRow(0),
  38   fSector(-1),
  39   fRelativeSector(-1),
  40   fCurrentSigmaY2(1e10),
  41   fCurrentSigmaZ2(1e10),
  42   fErrorY2(1e10),
  43   fErrorZ2(1e10),
  44   fCurrentCluster(0x0),
  45   fCurrentClusterIndex1(-1),
  46   fInDead(kFALSE),
  47   fIsSeeding(kFALSE),
  48   fNoCluster(0),
  49   fSort(0),
  50   fBSigned(kFALSE),
  51   fSeedType(0),
  52   fSeed1(-1),
  53   fSeed2(-1),
  54   fMAngular(0),
  55   fCircular(0)
  56 {
  57   //
  58   for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
  59   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
  60   for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
  61   for (Int_t i=0;i<AliPID::kSPECIES;i++)   fTPCr[i]=0.2;
  62   for (Int_t i=0;i<4;i++) {
  63     fDEDX[i] = 0.;
  64     fSDEDX[i] = 1e10;
  65     fNCDEDX[i] = 0;
  66   }
  67   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
  68 }
  69
  70 AliTPCseed::AliTPCseed(const AliTPCseed &s):
  71   AliTPCtrack(s),
  72   fEsd(0x0),
  73   fPoints(0x0),
  74   fEPoints(0x0),
  75   fRow(0),
  76   fSector(-1),
  77   fRelativeSector(-1),
  78   fCurrentSigmaY2(1e10),
  79   fCurrentSigmaZ2(1e10),
  80   fErrorY2(1e10),
  81   fErrorZ2(1e10),
  82   fCurrentCluster(0x0),
  83   fCurrentClusterIndex1(-1),
  84   fInDead(kFALSE),
  85   fIsSeeding(kFALSE),
  86   fNoCluster(0),
  87   fSort(0),
  88   fBSigned(kFALSE),
  89   fSeedType(0),
  90   fSeed1(-1),
  91   fSeed2(-1),
  92   fMAngular(0),
  93   fCircular(0)
  94 {
  95   //---------------------
  96   // dummy copy constructor
  97   //-------------------------
  98   for (Int_t i=0;i<160;i++) fClusterPointer[i] = s.fClusterPointer[i];
  99   for (Int_t i=0;i<160;i++) fIndex[i] = s.fIndex[i];
 100 }
 101 AliTPCseed::AliTPCseed(const AliTPCtrack &t):
 102   AliTPCtrack(t),
 103   fEsd(0x0),
 104   fPoints(0x0),
 105   fEPoints(0x0),
 106   fRow(0),
 107   fSector(-1),
 108   fRelativeSector(-1),
 109   fCurrentSigmaY2(1e10),
 110   fCurrentSigmaZ2(1e10),
 111   fErrorY2(1e10),
 112   fErrorZ2(1e10),
 113   fCurrentCluster(0x0),
 114   fCurrentClusterIndex1(-1),
 115   fInDead(kFALSE),
 116   fIsSeeding(kFALSE),
 117   fNoCluster(0),
 118   fSort(0),
 119   fBSigned(kFALSE),
 120   fSeedType(0),
 121   fSeed1(-1),
 122   fSeed2(-1),
 123   fMAngular(0),
 124   fCircular(0)
 125 {
 126   //
 127   // Constructor from AliTPCtrack
 128   //
 129   fFirstPoint =0;
 130   for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=t.GetKinkIndex(i);
 131   for (Int_t i=0;i<5;i++)   fTPCr[i]=0.2;
 132   for (Int_t i=0;i<160;i++) {
 133     fClusterPointer[i] = 0;
 134     Int_t index = t.GetClusterIndex(i);
 135     if (index>=-1){
 136       SetClusterIndex2(i,index);
 137     }
 138     else{
 139       SetClusterIndex2(i,-3);
 140     }
 141   }
 142   for (Int_t i=0;i<4;i++) {
 143     fDEDX[i] = 0.;
 144     fSDEDX[i] = 1e10;
 145     fNCDEDX[i] = 0;
 146   }
 147   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
 148 }
 149
 150 AliTPCseed::AliTPCseed(UInt_t index,  const Double_t xx[5],
 151                        const Double_t cc[15],
 152                        Double_t xr, Double_t alpha):
 153   AliTPCtrack(index, xx, cc, xr, alpha),
 154   fEsd(0x0),
 155   fPoints(0x0),
 156   fEPoints(0x0),
 157   fRow(0),
 158   fSector(-1),
 159   fRelativeSector(-1),
 160   fCurrentSigmaY2(1e10),
 161   fCurrentSigmaZ2(1e10),
 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 {
 177   //
 178   // Constructor
 179   //
 180   fFirstPoint =0;
 181   for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
 182   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
 183   for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
 184   for (Int_t i=0;i<5;i++)   fTPCr[i]=0.2;
 185   for (Int_t i=0;i<4;i++) {
 186     fDEDX[i] = 0.;
 187     fSDEDX[i] = 1e10;
 188     fNCDEDX[i] = 0;
 189   }
 190   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
 191 }
 192
 193 AliTPCseed::~AliTPCseed(){
 194   //
 195   // destructor
 196   if (fPoints) delete fPoints;
 197   fPoints =0;
 198   if (fEPoints) delete fEPoints;
 199   fEPoints = 0;
 200   fNoCluster =0;
 201 }
 202
 203 AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
 204 {
 205   //
 206   //
 207   return &fTrackPoints[i];
 208 }
 209
 210 void AliTPCseed::RebuildSeed()
 211 {
 212   //
 213   // rebuild seed to be ready for storing
 214   AliTPCclusterMI cldummy;
 215   cldummy.SetQ(0);
 216   AliTPCTrackPoint pdummy;
 217   pdummy.GetTPoint().fIsShared = 10;
 218   for (Int_t i=0;i<160;i++){
 219     AliTPCclusterMI * cl0 = fClusterPointer[i];
 220     AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
 221     if (cl0){
 222       trpoint->GetTPoint() = *(GetTrackPoint(i));
 223       trpoint->GetCPoint() = *cl0;
 224       trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
 225     }
 226     else{
 227       *trpoint = pdummy;
 228       trpoint->GetCPoint()= cldummy;
 229     }
 230
 231   }
 232
 233 }
 234
 235
 236 Double_t AliTPCseed::GetDensityFirst(Int_t n)
 237 {
 238   //
 239   //
 240   // return cluster for n rows bellow first point
 241   Int_t nfoundable = 1;
 242   Int_t nfound      = 1;
 243   for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
 244     Int_t index = GetClusterIndex2(i);
 245     if (index!=-1) nfoundable++;
 246     if (index>0) nfound++;
 247   }
 248   if (nfoundable<n) return 0;
 249   return Double_t(nfound)/Double_t(nfoundable);
 250
 251 }
 252
 253
 254 void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
 255 {
 256   // get cluster stat.  on given region
 257   //
 258   found       = 0;
 259   foundable   = 0;
 260   shared      =0;
 261   for (Int_t i=first;i<last; i++){
 262     Int_t index = GetClusterIndex2(i);
 263     if (index!=-1) foundable++;
 264     if (fClusterPointer[i]) {
 265       found++;
 266     }
 267     else
 268       continue;
 269
 270     if (fClusterPointer[i]->IsUsed(10)) {
 271       shared++;
 272       continue;
 273     }
 274     if (!plus2) continue; //take also neighborhoud
 275     //
 276     if ( (i>0) && fClusterPointer[i-1]){
 277       if (fClusterPointer[i-1]->IsUsed(10)) {
 278         shared++;
 279         continue;
 280       }
 281     }
 282     if ( fClusterPointer[i+1]){
 283       if (fClusterPointer[i+1]->IsUsed(10)) {
 284         shared++;
 285         continue;
 286       }
 287     }
 288
 289   }
 290   //if (shared>found){
 291     //Error("AliTPCseed::GetClusterStatistic","problem\n");
 292   //}
 293 }
 294
 295
 296
 297
 298
 299 void AliTPCseed::Reset(Bool_t all)
 300 {
 301   //
 302   //
 303   SetNumberOfClusters(0);
 304   fNFoundable = 0;
 305   SetChi2(0);
 306   ResetCovariance();
 307   /*
 308   if (fTrackPoints){
 309     for (Int_t i=0;i<8;i++){
 310       delete [] fTrackPoints[i];
 311     }
 312     delete fTrackPoints;
 313     fTrackPoints =0;
 314   }
 315   */
 316
 317   if (all){
 318     for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
 319     for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
 320   }
 321
 322 }
 323
 324
 325 void AliTPCseed::Modify(Double_t factor)
 326 {
 327
 328   //------------------------------------------------------------------
 329   //This function makes a track forget its history :)
 330   //------------------------------------------------------------------
 331   if (factor<=0) {
 332     ResetCovariance();
 333     return;
 334   }
 335   fC00*=factor;
 336   fC10*=0;  fC11*=factor;
 337   fC20*=0;  fC21*=0;  fC22*=factor;
 338   fC30*=0;  fC31*=0;  fC32*=0;  fC33*=factor;
 339   fC40*=0;  fC41*=0;  fC42*=0;  fC43*=0;  fC44*=factor;
 340   SetNumberOfClusters(0);
 341   fNFoundable =0;
 342   SetChi2(0);
 343   fRemoval = 0;
 344   fCurrentSigmaY2 = 0.000005;
 345   fCurrentSigmaZ2 = 0.000005;
 346   fNoCluster     = 0;
 347   //fFirstPoint = 160;
 348   //fLastPoint  = 0;
 349 }
 350
 351
 352
 353
 354 Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
 355 {
 356   //-----------------------------------------------------------------
 357   // This function find proloncation of a track to a reference plane x=xk.
 358   // doesn't change internal state of the track
 359   //-----------------------------------------------------------------
 360
 361   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
 362
 363   if (TMath::Abs(fP4*xk - fP2) >= AliTPCReconstructor::GetMaxSnpTrack()) {
 364     return 0;
 365   }
 366
 367   //  Double_t y1=fP0, z1=fP1;
 368   Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
 369   Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
 370
 371   y = fP0;
 372   z = fP1;
 373   //y += dx*(c1+c2)/(r1+r2);
 374   //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
 375
 376   Double_t dy = dx*(c1+c2)/(r1+r2);
 377   Double_t dz = 0;
 378   //
 379   Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
 380   /*
 381   if (TMath::Abs(delta)>0.0001){
 382     dz = fP3*TMath::ASin(delta)/fP4;
 383   }else{
 384     dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
 385   }
 386   */
 387   //  dz =  fP3*AliTPCFastMath::FastAsin(delta)/fP4;
 388   dz =  fP3*TMath::ASin(delta)/fP4;
 389   //
 390   y+=dy;
 391   z+=dz;
 392
 393
 394   return 1;
 395 }
 396
 397
 398 //_____________________________________________________________________________
 399 Double_t AliTPCseed::GetPredictedChi2(const AliCluster *c) const
 400 {
 401   //-----------------------------------------------------------------
 402   // This function calculates a predicted chi2 increment.
 403   //-----------------------------------------------------------------
 404   //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
 405   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
 406   r00+=fC00; r01+=fC10; r11+=fC11;
 407
 408   Double_t det=r00*r11 - r01*r01;
 409   if (TMath::Abs(det) < 1.e-10) {
 410     //Int_t n=GetNumberOfClusters();
 411     //if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
 412     return 1e10;
 413   }
 414   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
 415
 416   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
 417
 418   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
 419 }
 420
 421
 422 //_________________________________________________________________________________________
 423
 424
 425 Int_t AliTPCseed::Compare(const TObject *o) const {
 426   //-----------------------------------------------------------------
 427   // This function compares tracks according to the sector - for given sector according z
 428   //-----------------------------------------------------------------
 429   AliTPCseed *t=(AliTPCseed*)o;
 430
 431   if (fSort == 0){
 432     if (t->fRelativeSector>fRelativeSector) return -1;
 433     if (t->fRelativeSector<fRelativeSector) return 1;
 434     Double_t z2 = t->GetZ();
 435     Double_t z1 = GetZ();
 436     if (z2>z1) return 1;
 437     if (z2<z1) return -1;
 438     return 0;
 439   }
 440   else {
 441     Float_t f2 =1;
 442     f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
 443     if (t->fBConstrain) f2=1.2;
 444
 445     Float_t f1 =1;
 446     f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
 447
 448     if (fBConstrain)   f1=1.2;
 449
 450     if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
 451     else return +1;
 452   }
 453 }
 454
 455
 456
 457
 458 //_____________________________________________________________________________
 459 Int_t AliTPCseed::Update(const AliCluster *c, Double_t chisq, UInt_t /*index*/) {
 460   //-----------------------------------------------------------------
 461   // This function associates a cluster with this track.
 462   //-----------------------------------------------------------------
 463   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
 464
 465   r00+=fC00; r01+=fC10; r11+=fC11;
 466   Double_t det=r00*r11 - r01*r01;
 467   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 468
 469   Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
 470   Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
 471   Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
 472   Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
 473   Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
 474
 475   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
 476   Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
 477   if (TMath::Abs(cur*fX-eta) >= AliTPCReconstructor::GetMaxSnpTrack()) {
 478     return 0;
 479   }
 480
 481   fP0 += k00*dy + k01*dz;
 482   fP1 += k10*dy + k11*dz;
 483   fP2  = eta;
 484   fP3 += k30*dy + k31*dz;
 485   fP4  = cur;
 486
 487   Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
 488   Double_t c12=fC21, c13=fC31, c14=fC41;
 489
 490   fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
 491   fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
 492   fC40-=k00*c04+k01*c14;
 493
 494   fC11-=k10*c01+k11*fC11;
 495   fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
 496   fC41-=k10*c04+k11*c14;
 497
 498   fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
 499   fC42-=k20*c04+k21*c14;
 500
 501   fC33-=k30*c03+k31*c13;
 502   fC43-=k40*c03+k41*c13;
 503
 504   fC44-=k40*c04+k41*c14;
 505
 506   Int_t n=GetNumberOfClusters();
 507   //  fIndex[n]=index;
 508   SetNumberOfClusters(n+1);
 509   SetChi2(GetChi2()+chisq);
 510
 511   return 1;
 512 }
 513
 514
 515
 516 //_____________________________________________________________________________
 517 Float_t AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
 518   //-----------------------------------------------------------------
 519   // This funtion calculates dE/dX within the "low" and "up" cuts.
 520   //-----------------------------------------------------------------
 521
 522   Float_t amp[200];
 523   Float_t angular[200];
 524   Float_t weight[200];
 525   Int_t index[200];
 526   //Int_t nc = 0;
 527   //  TClonesArray & arr = *fPoints;
 528   Float_t meanlog = 100.;
 529
 530   Float_t mean[4]  = {0,0,0,0};
 531   Float_t sigma[4] = {1000,1000,1000,1000};
 532   Int_t nc[4]      = {0,0,0,0};
 533   Float_t norm[4]    = {1000,1000,1000,1000};
 534   //
 535   //
 536   fNShared =0;
 537
 538   for (Int_t of =0; of<4; of++){
 539     for (Int_t i=of+i1;i<i2;i+=4)
 540       {
 541         Int_t index = fIndex[i];
 542         if (index<0||index&0x8000) continue;
 543
 544         //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
 545         AliTPCTrackerPoint * point = GetTrackPoint(i);
 546         //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
 547         //AliTPCTrackerPoint * pointp = 0;
 548         //if (i<159) pointp = GetTrackPoint(i+1);
 549
 550         if (point==0) continue;
 551         AliTPCclusterMI * cl = fClusterPointer[i];
 552         if (cl==0) continue;
 553         if (onlyused && (!cl->IsUsed(10))) continue;
 554         if (cl->IsUsed(11)) {
 555           fNShared++;
 556           continue;
 557         }
 558         Int_t   type   = cl->GetType();
 559         //if (point->fIsShared){
 560         //  fNShared++;
 561         //  continue;
 562         //}
 563         //if (pointm)
 564         //  if (pointm->fIsShared) continue;
 565         //if (pointp)
 566         //  if (pointp->fIsShared) continue;
 567
 568         if (type<0) continue;
 569         //if (type>10) continue;
 570         //if (point->GetErrY()==0) continue;
 571         //if (point->GetErrZ()==0) continue;
 572
 573         //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
 574         //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
 575         //if ((ddy*ddy+ddz*ddz)>10) continue;
 576
 577
 578         //      if (point->GetCPoint().GetMax()<5) continue;
 579         if (cl->GetMax()<5) continue;
 580         Float_t angley = point->GetAngleY();
 581         Float_t anglez = point->GetAngleZ();
 582
 583         Float_t rsigmay2 =  point->GetSigmaY();
 584         Float_t rsigmaz2 =  point->GetSigmaZ();
 585         /*
 586         Float_t ns = 1.;
 587         if (pointm){
 588           rsigmay +=  pointm->GetTPoint().GetSigmaY();
 589           rsigmaz +=  pointm->GetTPoint().GetSigmaZ();
 590           ns+=1.;
 591         }
 592         if (pointp){
 593           rsigmay +=  pointp->GetTPoint().GetSigmaY();
 594           rsigmaz +=  pointp->GetTPoint().GetSigmaZ();
 595           ns+=1.;
 596         }
 597         rsigmay/=ns;
 598         rsigmaz/=ns;
 599         */
 600
 601         Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
 602
 603         Float_t ampc   = 0;     // normalization to the number of electrons
 604         if (i>64){
 605           //      ampc = 1.*point->GetCPoint().GetMax();
 606           ampc = 1.*cl->GetMax();
 607           //ampc = 1.*point->GetCPoint().GetQ();
 608           //      AliTPCClusterPoint & p = point->GetCPoint();
 609           //      Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
 610           // Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
 611           //Float_t dz =
 612           //  TMath::Abs( Int_t(iz) - iz + 0.5);
 613           //ampc *= 1.15*(1-0.3*dy);
 614           //ampc *= 1.15*(1-0.3*dz);
 615           //      Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
 616           //ampc               *=zfactor;
 617         }
 618         else{
 619           //ampc = 1.0*point->GetCPoint().GetMax();
 620           ampc = 1.0*cl->GetMax();
 621           //ampc = 1.0*point->GetCPoint().GetQ();
 622           //AliTPCClusterPoint & p = point->GetCPoint();
 623           // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
 624           //Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
 625           //Float_t dz =
 626           //  TMath::Abs( Int_t(iz) - iz + 0.5);
 627
 628           //ampc *= 1.15*(1-0.3*dy);
 629           //ampc *= 1.15*(1-0.3*dz);
 630           //    Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
 631           //ampc               *=zfactor;
 632
 633         }
 634         ampc *= 2.0;     // put mean value to channel 50
 635         //ampc *= 0.58;     // put mean value to channel 50
 636         Float_t w      =  1.;
 637         //      if (type>0)  w =  1./(type/2.-0.5);
 638         //      Float_t z = TMath::Abs(cl->GetZ());
 639         if (i<64) {
 640           ampc /= 0.6;
 641           //ampc /= (1+0.0008*z);
 642         } else
 643           if (i>128){
 644             ampc /=1.5;
 645             //ampc /= (1+0.0008*z);
 646           }else{
 647             //ampc /= (1+0.0008*z);
 648           }
 649
 650         if (type<0) {  //amp at the border - lower weight
 651           // w*= 2.;
 652
 653           continue;
 654         }
 655         if (rsigma>1.5) ampc/=1.3;  // if big backround
 656         amp[nc[of]]        = ampc;
 657         angular[nc[of]]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
 658         weight[nc[of]]     = w;
 659         nc[of]++;
 660       }
 661
 662     TMath::Sort(nc[of],amp,index,kFALSE);
 663     Float_t sumamp=0;
 664     Float_t sumamp2=0;
 665     Float_t sumw=0;
 666     //meanlog = amp[index[Int_t(nc[of]*0.33)]];
 667     meanlog = 50;
 668     for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
 669       Float_t ampl      = amp[index[i]]/angular[index[i]];
 670       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
 671       //
 672       sumw    += weight[index[i]];
 673       sumamp  += weight[index[i]]*ampl;
 674       sumamp2 += weight[index[i]]*ampl*ampl;
 675       norm[of]    += angular[index[i]]*weight[index[i]];
 676     }
 677     if (sumw<1){
 678       SetdEdx(0);
 679     }
 680     else {
 681       norm[of] /= sumw;
 682       mean[of]  = sumamp/sumw;
 683       sigma[of] = sumamp2/sumw-mean[of]*mean[of];
 684       if (sigma[of]>0.1)
 685         sigma[of] = TMath::Sqrt(sigma[of]);
 686       else
 687         sigma[of] = 1000;
 688
 689     mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
 690     //mean  *=(1-0.02*(sigma/(mean*0.17)-1.));
 691     //mean *=(1-0.1*(norm-1.));
 692     }
 693   }
 694
 695   Float_t dedx =0;
 696   fSdEdx =0;
 697   fMAngular =0;
 698   //  mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
 699   //  mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
 700
 701
 702   //  dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
 703   //  (  TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
 704
 705   Int_t norm2 = 0;
 706   Int_t norm3 = 0;
 707   for (Int_t i =0;i<4;i++){
 708     if (nc[i]>2&&nc[i]<1000){
 709       dedx      += mean[i] *nc[i];
 710       fSdEdx    += sigma[i]*(nc[i]-2);
 711       fMAngular += norm[i] *nc[i];
 712       norm2     += nc[i];
 713       norm3     += nc[i]-2;
 714     }
 715     fDEDX[i]  = mean[i];
 716     fSDEDX[i] = sigma[i];
 717     fNCDEDX[i]= nc[i];
 718   }
 719
 720   if (norm3>0){
 721     dedx   /=norm2;
 722     fSdEdx /=norm3;
 723     fMAngular/=norm2;
 724   }
 725   else{
 726     SetdEdx(0);
 727     return 0;
 728   }
 729   //  Float_t dedx1 =dedx;
 730   /*
 731   dedx =0;
 732   for (Int_t i =0;i<4;i++){
 733     if (nc[i]>2&&nc[i]<1000){
 734       mean[i]   = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
 735       dedx      += mean[i] *nc[i];
 736     }
 737     fDEDX[i]  = mean[i];
 738   }
 739   dedx /= norm2;
 740   */
 741
 742
 743   SetdEdx(dedx);
 744   return dedx;
 745 }
 746 Double_t AliTPCseed::Bethe(Double_t bg){
 747   //
 748   // This is the Bethe-Bloch function normalised to 1 at the minimum
 749   //
 750   Double_t bg2=bg*bg;
 751   Double_t bethe;
 752   if (bg<3.5e1)
 753     bethe=(1.+ bg2)/bg2*(log(5940*bg2) - bg2/(1.+ bg2));
 754   else // Density effect ( approximately :)
 755     bethe=1.15*(1.+ bg2)/bg2*(log(3.5*5940*bg) - bg2/(1.+ bg2));
 756   return bethe/11.091;
 757 }
 758
 759 void AliTPCseed::CookPID()
 760 {
 761   //
 762   // cook PID information according dEdx
 763   //
 764   Double_t fRange = 10.;
 765   Double_t fRes   = 0.1;
 766   Double_t fMIP   = 47.;
 767   //
 768   Int_t ns=AliPID::kSPECIES;
 769   Double_t sumr =0;
 770   for (Int_t j=0; j<ns; j++) {
 771     Double_t mass=AliPID::ParticleMass(j);
 772     Double_t mom=P();
 773     Double_t dedx=fdEdx/fMIP;
 774     Double_t bethe=Bethe(mom/mass);
 775     Double_t sigma=fRes*bethe;
 776     if (sigma>0.001){
 777       if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 778         fTPCr[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 779         sumr+=fTPCr[j];
 780         continue;
 781       }
 782       fTPCr[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 783       sumr+=fTPCr[j];
 784     }
 785     else{
 786       fTPCr[j]=1.;
 787       sumr+=fTPCr[j];
 788     }
 789   }
 790   for (Int_t j=0; j<ns; j++) {
 791     fTPCr[j]/=sumr;           //normalize
 792   }
 793 }
 794
 795 /*
 796 void AliTPCseed::CookdEdx2(Double_t low, Double_t up) {
 797   //-----------------------------------------------------------------
 798   // This funtion calculates dE/dX within the "low" and "up" cuts.
 799   //-----------------------------------------------------------------
 800
 801   Float_t amp[200];
 802   Float_t angular[200];
 803   Float_t weight[200];
 804   Int_t index[200];
 805   Bool_t inlimit[200];
 806   for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
 807   for (Int_t i=0;i<200;i++) amp[i]=10000;
 808   for (Int_t i=0;i<200;i++) angular[i]= 1;;
 809
 810
 811   //
 812   Float_t meanlog = 100.;
 813   Int_t indexde[4]={0,64,128,160};
 814
 815   Float_t amean     =0;
 816   Float_t asigma    =0;
 817   Float_t anc       =0;
 818   Float_t anorm     =0;
 819
 820   Float_t mean[4]  = {0,0,0,0};
 821   Float_t sigma[4] = {1000,1000,1000,1000};
 822   Int_t nc[4]      = {0,0,0,0};
 823   Float_t norm[4]    = {1000,1000,1000,1000};
 824   //
 825   //
 826   fNShared =0;
 827
 828   //  for (Int_t of =0; of<3; of++){
 829   //  for (Int_t i=indexde[of];i<indexde[of+1];i++)
 830   for (Int_t i =0; i<160;i++)
 831     {
 832         AliTPCTrackPoint * point = GetTrackPoint(i);
 833         if (point==0) continue;
 834         if (point->fIsShared){
 835           fNShared++;
 836           continue;
 837         }
 838         Int_t   type   = point->GetCPoint().GetType();
 839         if (type<0) continue;
 840         if (point->GetCPoint().GetMax()<5) continue;
 841         Float_t angley = point->GetTPoint().GetAngleY();
 842         Float_t anglez = point->GetTPoint().GetAngleZ();
 843         Float_t rsigmay =  point->GetCPoint().GetSigmaY();
 844         Float_t rsigmaz =  point->GetCPoint().GetSigmaZ();
 845         Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
 846
 847         Float_t ampc   = 0;     // normalization to the number of electrons
 848         if (i>64){
 849           ampc =  point->GetCPoint().GetMax();
 850         }
 851         else{
 852           ampc = point->GetCPoint().GetMax();
 853         }
 854         ampc *= 2.0;     // put mean value to channel 50
 855         //      ampc *= 0.565;     // put mean value to channel 50
 856
 857         Float_t w      =  1.;
 858         Float_t z = TMath::Abs(point->GetCPoint().GetZ());
 859         if (i<64) {
 860           ampc /= 0.63;
 861         } else
 862           if (i>128){
 863             ampc /=1.51;
 864           }
 865         if (type<0) {  //amp at the border - lower weight
 866           continue;
 867         }
 868         if (rsigma>1.5) ampc/=1.3;  // if big backround
 869         angular[i]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
 870         amp[i]        = ampc/angular[i];
 871         weight[i]     = w;
 872         anc++;
 873     }
 874
 875   TMath::Sort(159,amp,index,kFALSE);
 876   for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){
 877     inlimit[index[i]] = kTRUE;  // take all clusters
 878   }
 879
 880   //  meanlog = amp[index[Int_t(anc*0.3)]];
 881   meanlog =10000.;
 882   for (Int_t of =0; of<3; of++){
 883     Float_t sumamp=0;
 884     Float_t sumamp2=0;
 885     Float_t sumw=0;
 886    for (Int_t i=indexde[of];i<indexde[of+1];i++)
 887       {
 888         if (inlimit[i]==kFALSE) continue;
 889         Float_t ampl      = amp[i];
 890         ///angular[i];
 891         ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
 892         //
 893         sumw    += weight[i];
 894         sumamp  += weight[i]*ampl;
 895         sumamp2 += weight[i]*ampl*ampl;
 896         norm[of]    += angular[i]*weight[i];
 897         nc[of]++;
 898       }
 899    if (sumw<1){
 900      SetdEdx(0);
 901    }
 902    else {
 903      norm[of] /= sumw;
 904      mean[of]  = sumamp/sumw;
 905      sigma[of] = sumamp2/sumw-mean[of]*mean[of];
 906      if (sigma[of]>0.1)
 907        sigma[of] = TMath::Sqrt(sigma[of]);
 908      else
 909        sigma[of] = 1000;
 910      mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
 911    }
 912   }
 913
 914   Float_t dedx =0;
 915   fSdEdx =0;
 916   fMAngular =0;
 917   //
 918   Int_t norm2 = 0;
 919   Int_t norm3 = 0;
 920   Float_t www[3] = {12.,14.,17.};
 921   //Float_t www[3] = {1.,1.,1.};
 922
 923   for (Int_t i =0;i<3;i++){
 924     if (nc[i]>2&&nc[i]<1000){
 925       dedx      += mean[i] *nc[i]*www[i]/sigma[i];
 926       fSdEdx    += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
 927       fMAngular += norm[i] *nc[i];
 928       norm2     += nc[i]*www[i]/sigma[i];
 929       norm3     += (nc[i]-2)*www[i]/sigma[i];
 930     }
 931     fDEDX[i]  = mean[i];
 932     fSDEDX[i] = sigma[i];
 933     fNCDEDX[i]= nc[i];
 934   }
 935
 936   if (norm3>0){
 937     dedx   /=norm2;
 938     fSdEdx /=norm3;
 939     fMAngular/=norm2;
 940   }
 941   else{
 942     SetdEdx(0);
 943     return;
 944   }
 945   //  Float_t dedx1 =dedx;
 946
 947   dedx =0;
 948   Float_t norm4 = 0;
 949   for (Int_t i =0;i<3;i++){
 950     if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
 951       //mean[i]   = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
 952       dedx      += mean[i] *(nc[i])/(sigma[i]);
 953       norm4     += (nc[i])/(sigma[i]);
 954     }
 955     fDEDX[i]  = mean[i];
 956   }
 957   if (norm4>0) dedx /= norm4;
 958
 959
 960
 961   SetdEdx(dedx);
 962
 963   //mi deDX
 964
 965 }
 966 */