TPC/AliTPCTrackHits.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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Time Projection Chamber  track hits object                                //
  21 //
  22 //  Origin: Marian Ivanov , GSI Darmstadt
  23 //
  24 //  Class for storing simulated AliTPCHits  for given track                  //
  25 //             -average compression comparing to classical ClonesArray is    //
  26 //              around 5-7 (depending on the required hit precision)       //
  27 //                                                                           //
  28 //Begin_Html
  29 /*
  30 <img src="gif/AliTPCTrackHits.gif">
  31 */
  32 //End_Html
  33 //                                                                           //
  34 //                                                                          //
  35 ///////////////////////////////////////////////////////////////////////////////
  36
  37 #include <TError.h>
  38
  39 #include "AliTPC.h"
  40 #include "AliTPCTrackHits.h"
  41 //           Interface classes
  42 #include "AliTPCTrackHitsInterfaces.h"
  43
  44 ClassImp(AliTPCTrackHits)
  45 LClassImp(AliTrackHitsInfo)
  46 LClassImp(AliTrackHitsParam)
  47 LClassImp(AliHitInfo)
  48
  49 Int_t AliTrackHitsInfo::fgCounter1 =0;
  50 Int_t AliTrackHitsInfo::fgCounter2 =0;
  51 Int_t AliTrackHitsParam::fgCounter1 =0;
  52 Int_t AliTrackHitsParam::fgCounter2 =0;
  53 Int_t AliHitInfo::fgCounter1 =0;
  54 Int_t AliHitInfo::fgCounter2 =0;
  55 Int_t AliTPCTrackHits::fgCounter1 =0;
  56 Int_t AliTPCTrackHits::fgCounter2 =0;
  57 const Double_t AliTPCTrackHits::fgkPrecision=1e-6;  //precision
  58 const Double_t AliTPCTrackHits::fgkPrecision2=1e-20;  //precision
  59 const Double_t AliTPCTrackHits::fgkTimePrecision=20.e-9;  //hit time precision
  60
  61
  62 class AliTPCCurrentHit {
  63 private:
  64   AliTPChit fHit;     //   - hit in "standard" representation
  65   UInt_t   fInfoIndex;//   - current info pointer
  66   UInt_t   fParamIndex;//  - current param pointer
  67   UInt_t   fStackIndex; // - current hit stack index
  68   Double_t fR;   //current Radius
  69   Bool_t  fStatus; //current status
  70 };
  71
  72
  73 class  AliTPCTempHitInfo {
  74 private:
  75   enum    { kStackSize = 100};
  76   AliTPCTempHitInfo();
  77   AliTPCTempHitInfo(const AliTPCTempHitInfo &)
  78     {::Fatal("copy ctor","Not implemented\n");}
  79   AliTPCTempHitInfo & operator = (const AliTPCTempHitInfo &)
  80     {::Fatal("= operator","Not implemented\n");return *this;}
  81
  82   void     NewParam(Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time);
  83   void     SetHit(Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time);
  84   Double_t * GetPosition(Int_t index){return &fPositionStack[index*3];}
  85   void    UpdateParam(Double_t maxdelta); //recal
  86   void   Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
  87             Double_t fSumX,  Double_t fSumX2, Double_t fSumX3,
  88             Double_t fSumX4, Int_t n,
  89               Double_t &a, Double_t &b, Double_t &c);
  90   void  Fit(AliTrackHitsParam * param);
  91   Double_t fSumDr;    //fSumDr
  92   Double_t fSumDr2;   //fSumDr2
  93   Double_t fSumDr3;   //  fSumDr3
  94   Double_t fSumDr4;   //fSumDr4
  95   Double_t fSumDFi;  //fSumDFi
  96   Double_t fSumDFiDr; // fSumDFiDr
  97   Double_t fSumDFiDr2;//fSumDFiDr2
  98   Double_t fSumDZ;     //fSumDZ
  99   Double_t fSumDZDr;  //fSumDZDr
 100   Double_t fSumDZDr2;  //fSumDZDr2
 101   Double_t fOldR;     //previos r
 102   Double_t fPositionStack[3*kStackSize];  //position stack
 103   UInt_t   fQStack[kStackSize];           //Q stack
 104   Float_t  fTimeStack[kStackSize];        //time stack
 105   UInt_t fStackIndex;   //current stack index
 106   UInt_t fInfoIndex;    //current track info index
 107   UInt_t fParamIndex;   //current track parameters index
 108   AliTrackHitsInfo  * fInfo; //current track info
 109   AliTrackHitsParam * fParam; //current track param
 110 };
 111
 112
 113 AliTPCTempHitInfo::AliTPCTempHitInfo()
 114 {
 115   //
 116   //set to default value
 117   fSumDr=fSumDr2=fSumDr3=fSumDr4=
 118     fSumDFi=fSumDFiDr=fSumDFiDr2=
 119     fSumDZ=fSumDZDr=fSumDZDr2=0;
 120   fStackIndex = 0;
 121   fInfoIndex  = 0;
 122   fParamIndex = 0;
 123 }
 124
 125
 126 void AliTPCTempHitInfo::NewParam(Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time)
 127 {
 128   //
 129   //reset stack and sum parameters
 130   //store line initial point
 131   fSumDr=fSumDr2=fSumDr3=fSumDr4=
 132     fSumDFi=fSumDFiDr=fSumDFiDr2=
 133     fSumDZ=fSumDZDr=fSumDZDr2=0;
 134   fStackIndex=0;
 135   fParam->fR = r;
 136   fOldR = r;
 137   fParam->fZ = z;
 138   fParam->fFi = fi;
 139   fParam->fAn = 0.;
 140   fParam->fAd = 0.;
 141   fParam->fTheta =0.;
 142   fParam->fThetaD =0.;
 143   SetHit(r,z,fi,q,time);
 144 }
 145
 146 void AliTPCTempHitInfo::SetHit(Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time)
 147 {
 148   //
 149   //add hit to the stack
 150   //recalculate new estimete of line parameters
 151   Double_t *f = GetPosition(fStackIndex);
 152   f[0] = r;
 153   f[1] = z;
 154   f[2] = fi;
 155   fQStack[fStackIndex]=q;
 156   fTimeStack[fStackIndex]=time;
 157   if (fStackIndex==0) return;
 158   Double_t dr  = (r-fParam->fR);
 159   if (TMath::Abs(dr)<AliTPCTrackHits::fgkPrecision) dr =AliTPCTrackHits::fgkPrecision;
 160   Double_t dfi = fi-fParam->fFi;
 161   Double_t dz  = z -fParam->fZ;
 162   Double_t dr2 =dr*dr;
 163   Double_t dr3 =dr2*dr;
 164   Double_t dr4 =dr3*dr;
 165   fSumDr +=dr;
 166   fSumDr2+=dr2;
 167   fSumDr3+=dr3;
 168   fSumDr4+=dr4;
 169   fSumDFi +=dfi;
 170   fSumDFiDr+=dfi*dr;
 171   fSumDFiDr2+=dfi*dr2;
 172   fSumDZ +=dz;
 173   fSumDZDr+=dz*dr;
 174   fSumDZDr2+=dz*dr2;
 175
 176   //update fit parameters
 177   //
 178   Double_t det = fSumDr2*fSumDr4-fSumDr3*fSumDr3;
 179   if (TMath::Abs(det)<AliTPCTrackHits::fgkPrecision2) return;
 180   if ( ( fStackIndex>1 )  ){
 181     fParam->fAn = (fSumDr4*fSumDFiDr-fSumDr3*fSumDFiDr2)/det;
 182     fParam->fAd = (fSumDr2*fSumDFiDr2-fSumDr3*fSumDFiDr)/det;
 183   }
 184   else
 185     fParam->fAn = fSumDFiDr/fSumDr2;
 186   if ( ( fStackIndex>1 )  ){
 187     fParam->fTheta = (fSumDr4*fSumDZDr-fSumDr3*fSumDZDr2)/det;
 188     fParam->fThetaD= (fSumDr2*fSumDZDr2-fSumDr3*fSumDZDr)/det;
 189   }
 190   else
 191     fParam->fTheta = fSumDZDr/fSumDr2;
 192 }
 193
 194
 195 void   AliTPCTempHitInfo::UpdateParam(Double_t maxdelta)
 196 {
 197   //recalc parameters not fixing origin point
 198   if (fStackIndex>5){
 199     Double_t a,b,c;
 200     a=b=c=0;
 201     Fit2(fSumDFi, fSumDFiDr, fSumDFiDr2, fSumDr,fSumDr2,fSumDr3,fSumDr4,
 202          fStackIndex, a,b,c);
 203     if (TMath::Abs(a)<maxdelta){
 204       fParam->fFi +=a/fParam->fR;
 205       fParam->fAn = b;
 206       fParam->fAd = c;
 207     }
 208     Fit2(fSumDZ, fSumDZDr, fSumDZDr2, fSumDr,fSumDr2,fSumDr3,fSumDr4,
 209          fStackIndex, a,b,c) ;
 210     if (TMath::Abs(a)<maxdelta){
 211       fParam->fZ +=a;
 212       fParam->fTheta = b;
 213       fParam->fThetaD = c;
 214     }
 215   }
 216
 217 }
 218 void   AliTPCTempHitInfo::Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
 219             Double_t fSumX,  Double_t fSumX2, Double_t fSumX3,
 220             Double_t fSumX4, Int_t n,
 221             Double_t &a, Double_t &b, Double_t &c)
 222 {
 223   //fit of second order
 224   Double_t det =
 225     n* (fSumX2*fSumX4-fSumX3*fSumX3) -
 226     fSumX*      (fSumX*fSumX4-fSumX3*fSumX2)+
 227     fSumX2*     (fSumX*fSumX3-fSumX2*fSumX2);
 228
 229   if (TMath::Abs(det)> AliTPCTrackHits::fgkPrecision) {
 230     a =
 231       (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
 232        fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
 233        fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det;
 234     b=
 235       (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
 236       fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
 237       fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
 238     c=
 239       (n*(fSumX2*fSumYX2-fSumYX*fSumX3)-
 240        fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
 241        fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;
 242   }
 243 }
 244
 245 void   AliTPCTempHitInfo::Fit(AliTrackHitsParam * param)
 246 {
 247   // fit fixing first and the last point
 248   //result stored in new param
 249   Double_t dx2  = (GetPosition(fStackIndex))[0]-fParam->fR;
 250   Double_t det = fSumDr4+dx2*fSumDr2-2*dx2*fSumDr3;
 251   if ( (TMath::Abs(det)> AliTPCTrackHits::fgkPrecision) &&
 252        ((TMath::Abs(dx2)> AliTPCTrackHits::fgkPrecision))){
 253     Double_t dfi2 = (GetPosition(fStackIndex))[1]-fParam->fFi;
 254     param->fAd = (fSumDFiDr2+dfi2*fSumDr-dx2*fSumDFiDr-dfi2*fSumDr3/dx2)/det;
 255     param->fAn  = (dfi2-param->fAd*dx2*dx2)/dx2;
 256
 257     Double_t dz2 = (GetPosition(fStackIndex))[1]-fParam->fZ;
 258     param->fTheta = (fSumDZDr2+dz2*fSumDr-dx2*fSumDZDr-dz2*fSumDr3/dx2)/det;
 259     param->fTheta  = (dz2-param->fAd*dx2*dx2)/dx2;
 260   }
 261
 262 }
 263
 264 //______________________________________________________________________
 265 AliTrackHitsInfo::AliTrackHitsInfo() :
 266   fTrackID(0),
 267   fVolumeID(0),
 268   fHitParamIndex(0)
 269 {
 270   //
 271   // Default constructor
 272   //
 273   fgCounter1++;
 274   fgCounter2++;
 275 }
 276
 277 //______________________________________________________________________
 278 AliTrackHitsParam::AliTrackHitsParam() :
 279   fR(0),
 280   fZ(0),
 281   fFi(0),
 282   fAn(0),
 283   fAd(0),
 284   fTheta(0),
 285   fThetaD(0)
 286 {
 287   //
 288   // Default constructor
 289   //
 290   fgCounter1++;
 291   fgCounter2++;
 292 }
 293
 294
 295 AliTPCTrackHits::AliTPCTrackHits()
 296 {
 297   //
 298   //default constructor
 299   //
 300   const Float_t kHitPrecision=0.002; //default precision for hit position in cm
 301   const Float_t kStep =0.003;  //30 mum step
 302   const UShort_t kMaxDistance =100;  //maximum distance 100
 303
 304   fPrecision=kHitPrecision; //precision in cm
 305   fStep = kStep; //step size
 306   fMaxDistance = kMaxDistance; //maximum distance
 307   fTempInfo =0;
 308   fTrackHitsInfo = new AliObjectArray("AliTrackHitsInfo");
 309   fTrackHitsParam = new AliObjectArray("AliTrackHitsParam");
 310   fHitsPosAndQ = new TArrayOfArrayVStack("AliHitInfo");
 311   fCurrentHit = new AliTPCCurrentHit;
 312   fgCounter1++;
 313   fgCounter2++;
 314
 315 }
 316
 317 AliTPCTrackHits::AliTPCTrackHits(const AliTPCTrackHits& r) : TObject(r)
 318 {
 319   //dummy
 320 }
 321 AliTPCTrackHits &AliTPCTrackHits::operator=(const AliTPCTrackHits& /* r */)
 322 {
 323   //dummy
 324   return *this;
 325 }
 326
 327 AliTPCTrackHits::~AliTPCTrackHits()
 328 {
 329   //
 330   //default destructor
 331   //
 332   if (fTrackHitsInfo) delete fTrackHitsInfo;
 333   if (fTrackHitsParam) delete fTrackHitsParam;
 334   if (fHitsPosAndQ) delete fHitsPosAndQ;
 335   if (fCurrentHit) delete fCurrentHit;
 336   if (fTempInfo) delete fTempInfo;
 337   fgCounter1--;
 338 }
 339
 340 void AliTPCTrackHits::Clear()
 341 {
 342   //
 343   //clear object
 344   fTrackHitsInfo->Clear();
 345   fTrackHitsParam->Clear();
 346   //fTrackHitsInfo->Resize(0);
 347   //fTrackHitsParam->Resize(0);
 348   fHitsPosAndQ->Clear();
 349
 350   if (fTempInfo){
 351     delete fTempInfo;
 352     fTempInfo =0;
 353   }
 354 }
 355
 356
 357 void AliTPCTrackHits::AddHitKartez(Int_t volumeID, Int_t trackID, Double_t x,
 358               Double_t y, Double_t z,Int_t q, Float_t time)
 359 {
 360   //add hits (cartesian)
 361   Double_t r = TMath::Sqrt(x*x+y*y);
 362   Double_t fi = TMath::ACos(x/r);
 363   if (y<0) fi*=-1.;
 364     AddHit(volumeID,trackID,r,z,fi,q,time);
 365 }
 366
 367 void AliTPCTrackHits::AddHit(Int_t volumeID, Int_t trackID,
 368                              Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time)
 369 {
 370   //
 371   Bool_t diff=kFALSE;
 372   if (!fTempInfo) { //initialsation of track
 373     fTempInfo = new AliTPCTempHitInfo;
 374     //
 375     if (fTrackHitsInfo->GetCapacity()<10) fTrackHitsInfo->Reserve(10);
 376     fTrackHitsInfo->Resize(1);
 377     fTempInfo->fInfoIndex =0;
 378     if (fTrackHitsParam->GetCapacity()<100) fTrackHitsParam->Reserve(100);
 379     fTrackHitsParam->Resize(1);
 380     //
 381     fTempInfo->fInfo =
 382       (AliTrackHitsInfo*) (fTrackHitsInfo->At(0));
 383     fTempInfo->fInfo->fVolumeID = volumeID;
 384     fTempInfo->fInfo->fTrackID = trackID;
 385     fTempInfo->fInfo->fHitParamIndex =0;
 386     fTempInfo->fInfoIndex = 0;
 387     //
 388     fTempInfo->fParam =
 389       (AliTrackHitsParam*) (fTrackHitsParam->At(0));
 390     fTempInfo->fParamIndex = 0;
 391     fTempInfo->NewParam(r,z,fi,q,time);
 392     return;
 393   }
 394
 395   Int_t size = fHitsPosAndQ->GetSize();
 396   if (size>(Int_t)fTempInfo->fParamIndex) {
 397     fTempInfo->fParamIndex++;
 398     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize())
 399       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);
 400     fTempInfo->fParam =
 401       (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));
 402     fTempInfo->NewParam(r,z,fi,q,time);
 403     return;
 404   }
 405
 406
 407   // if new volume or new trackID
 408   if ( (volumeID!=fTempInfo->fInfo->fVolumeID) ||
 409        (trackID!=fTempInfo->fInfo->fTrackID)){
 410     diff=kTRUE;
 411
 412     FlushHitStack(kTRUE);
 413
 414     fTempInfo->fInfoIndex++;
 415     if (fTempInfo->fInfoIndex+1>fTrackHitsInfo->GetSize())
 416       fTrackHitsInfo->Resize(fTempInfo->fInfoIndex+1);
 417     fTempInfo->fInfo =
 418       (AliTrackHitsInfo*) (fTrackHitsInfo->At(fTempInfo->fInfoIndex));
 419     fTempInfo->fInfo->fVolumeID = volumeID;
 420     fTempInfo->fInfo->fTrackID = trackID;
 421     fTempInfo->fInfo->fHitParamIndex =fTempInfo->fParamIndex+1;
 422     // FlushHitStack(kTRUE);
 423
 424     fTempInfo->fParamIndex++;
 425     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize())
 426       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);
 427     fTempInfo->fParam =
 428       (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));
 429     fTempInfo->NewParam(r,z,fi,q,time);
 430     return;
 431   }
 432
 433   //calculate current fit precission to next point
 434   AliTrackHitsParam &param = *(fTempInfo->fParam);
 435   Double_t dd=0;
 436   Double_t dl=0;
 437   Double_t ratio=0;
 438   Double_t dr,dz,dfi,ddz,ddfi;
 439   Double_t drhit,ddl;
 440   dr=dz=dfi=ddz=ddfi=0;
 441   drhit = r-fTempInfo->fOldR;
 442   {
 443     //Double_t dfi2 = param.fAn+2*param.fAd*(r-param.fR);
 444     Double_t dfi2 = param.fAn;
 445     dfi2*=dfi2*fTempInfo->fOldR*fTempInfo->fOldR;
 446     //Double_t ddz2 =  param.fTheta+2*param.fThetaD*(r-param.fR);
 447     Double_t ddz2 =  param.fTheta;
 448     ddz2*=ddz2;
 449     ratio = TMath::Sqrt(1.+ dfi2+ ddz2);
 450   }
 451   //
 452
 453
 454   dl = (TMath::Abs(drhit*ratio/fStep)<32000) ?  fStep * Short_t(TMath::Nint(drhit*ratio/fStep)):0;
 455   ddl = dl - drhit*ratio;
 456   fTempInfo->fOldR += dl/ratio;
 457
 458   if (fTempInfo->fStackIndex>2){
 459     dr = r-param.fR;
 460     dz =  z-param.fZ;
 461     dfi = fi-param.fFi;
 462     ddz = dr*param.fTheta+dr*dr*param.fThetaD-dz;
 463     ddfi= dr*param.fAn+dr*dr*param.fAd-dfi;
 464     dd  = TMath::Sqrt(ddz*ddz+r*r*ddfi*ddfi+ddl*ddl);
 465     //
 466   }
 467   //safety factor 1.25
 468   if ( ( (dd*1.25>fPrecision) ) ||
 469        (fTempInfo->fStackIndex+4>fTempInfo->kStackSize) ||
 470        (TMath::Abs(dl/fStep)>fMaxDistance)  )
 471     diff=kTRUE;
 472   else{
 473     fTempInfo->fStackIndex++;
 474     fTempInfo->SetHit(r,z,fi,q,time);
 475     return;
 476   }
 477   //if parameter changed
 478   if (FlushHitStack(kFALSE)){   //if full buffer flushed
 479     fTempInfo->fParamIndex++;
 480     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize())
 481       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);
 482     fTempInfo->fParam =
 483       (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));
 484     fTempInfo->NewParam(r,z,fi,q,time);
 485   }
 486   else{
 487     fTempInfo->fStackIndex++;
 488     fTempInfo->SetHit(r,z,fi,q,time);
 489   }
 490 }
 491
 492 Bool_t AliTPCTrackHits::FlushHitStack(Bool_t force)
 493 {
 494   //
 495   //write fHitsPosAndQ information from the stack to te arrays
 496   if (!fTempInfo) return kFALSE;
 497   Int_t size = fHitsPosAndQ->GetSize();
 498
 499   if ( (size>0)&&(size!=(Int_t)fTempInfo->fParamIndex)) return kFALSE;
 500
 501   if (fHitsPosAndQ->Push(fTempInfo->fStackIndex+1)!=fTempInfo->fParamIndex){
 502     cout<<"internal error - contact MI\n";
 503     return kFALSE;
 504   }
 505   AliHitInfo * info;
 506
 507   AliTrackHitsParam & param = *(fTempInfo->fParam);
 508   //recalculate track parameter not fixing first point
 509   fTempInfo->UpdateParam(fStep/4.);
 510   //fTempInfo->Fit(fTempInfo->fParam);  //- fixing the first and the last point
 511
 512   Double_t oldr = param.fR;
 513   //cout<<"C3"<<fTempInfo->fStackIndex<<"\n"<<flush;
 514   UInt_t i;
 515   Double_t dd;
 516   for (i=0; i <= fTempInfo->fStackIndex; i++){
 517     Double_t * position = fTempInfo->GetPosition(i);
 518     Double_t   dr = position[0]-oldr;
 519     Double_t   ratio;
 520     {
 521       //Double_t dfi2 = param.fAn+2*param.fAd*(position[0]-param.fR);
 522       Double_t dfi2 = param.fAn;
 523       dfi2*=dfi2*oldr*oldr;
 524       //Double_t ddz2 =  param.fTheta+2*param.fThetaD*(position[0]-param.fR);
 525       Double_t ddz2 =  param.fTheta;
 526       ddz2*=ddz2;
 527       ratio = TMath::Sqrt(1.+ dfi2+ ddz2);
 528     }
 529
 530     Double_t dl = (TMath::Abs(dr*ratio/fStep)<32000) ? fStep*(Short_t)TMath::Nint(dr*ratio/fStep):0;
 531     dr = dl/ratio;
 532     oldr+=dr;
 533     //calculate precission
 534     AliTrackHitsParam &param = *(fTempInfo->fParam);
 535     //real deltas
 536     Double_t dr1=  position[0]-param.fR;
 537     Double_t dz =  position[1]-param.fZ;
 538
 539     Double_t dfi = position[2]-param.fFi;
 540     //extrapolated deltas
 541     Double_t dr2 = oldr-param.fR;
 542     Double_t ddr = dr2-dr1;
 543     Double_t ddz = dr2*param.fTheta+dr2*dr2*param.fThetaD-dz;
 544     Double_t ddfi= dr2*param.fAn+dr2*dr2*param.fAd-dfi;
 545     dd = TMath::Sqrt(ddz*ddz+oldr*oldr*ddfi*ddfi+ddr*ddr);
 546
 547     if ( (dd>fPrecision) ){
 548       if (i==0){
 549         param.fAn = 0;
 550         param.fAd = 0;
 551         param.fTheta =0;
 552         param.fThetaD =0;
 553         Double_t ddz = dr2*param.fTheta+dr2*dr2*param.fThetaD-dz;
 554         Double_t ddfi= dr2*param.fAn+dr2*dr2*param.fAd-dfi;
 555         dl = 0;
 556         dd = TMath::Sqrt(ddz*ddz+oldr*oldr*ddfi*ddfi+ddr*ddr);
 557       }
 558       else
 559         break;
 560     }
 561
 562     info = (AliHitInfo*)(fHitsPosAndQ->At(fTempInfo->fParamIndex,i));
 563     info->fHitDistance = (TMath::Abs(dl/fStep)<32000) ?Short_t(TMath::Nint(dl/fStep)):0;
 564     info->fCharge = Short_t(fTempInfo->fQStack[i]);
 565     info->fTime = TMath::Nint(fTempInfo->fTimeStack[i]/AliTPCTrackHits::fgkTimePrecision);
 566     /*
 567     cout<<"C2";
 568     cout<<" "<<fTempInfo->fStackIndex<<" \t";
 569     cout<<" "<<i<<" \t";
 570     cout<<position[0]<<"\t";
 571     cout<<position[1]<<"\t";
 572     cout<<position[2]<<"\t";
 573     cout<<param.fAn<<"\t";
 574     cout<<param.fTheta<<"\t";
 575     cout<<dr1<<"\t"<<ddr<<"\t"<<ddz<<"\t"<<ddfi<<"\t"<<dd<<"\n"<<flush;
 576     */
 577   }
 578
 579   if (i<=fTempInfo->fStackIndex){ //if previous iteration not succesfull
 580     fHitsPosAndQ->Resize(fTempInfo->fParamIndex,i);
 581     //
 582     fTempInfo->fParamIndex++;
 583     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize())
 584       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);
 585     fTempInfo->fParam =
 586         (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));
 587     Double_t * p = fTempInfo->GetPosition(i);
 588     UInt_t index2 = fTempInfo->fStackIndex;
 589     fTempInfo->NewParam(p[0],p[1],p[2],fTempInfo->fQStack[i],fTempInfo->fTimeStack[i]);
 590     if (i+1<=index2) FlushHitStack2(i+1,index2);
 591
 592     if (force) return      FlushHitStack(kTRUE);
 593     return kFALSE;
 594   }
 595   return kTRUE;
 596 }
 597
 598
 599 void AliTPCTrackHits::FlushHitStack2(Int_t index1, Int_t index2)
 600 {
 601   //
 602   // second iteration flush stack
 603   // call only for hits where first iteration were not succesfully interpolated
 604   Double_t * positionstack = new Double_t[3*(index2-index1+1)];
 605   UInt_t   * qstack        = new UInt_t[index2-index1+1];
 606   Float_t  * timestack        = new Float_t[index2-index1+1];
 607   memcpy(positionstack, &fTempInfo->fPositionStack[3*index1],
 608          (3*(index2-index1+1))*sizeof(Double_t));
 609   memcpy(qstack, &fTempInfo->fQStack[index1],(index2-index1+1)*sizeof(UInt_t));
 610   memcpy(timestack, &fTempInfo->fTimeStack[index1],(index2-index1+1)*sizeof(Float_t));
 611   Double_t *p = positionstack;
 612   for (Int_t j=0; j<=index2-index1;j++){
 613     fTempInfo->fStackIndex++;
 614     fTempInfo->SetHit(p[3*j+0],p[3*j+1],p[3*j+2],qstack[j],timestack[j]);
 615   }
 616   delete []positionstack;
 617   delete []qstack;
 618   delete []timestack;
 619 }
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629 Bool_t AliTPCTrackHits::First()
 630 {
 631   //
 632   //set Current hit for the first hit
 633   //
 634   AliTrackHitsInfo *info = (AliTrackHitsInfo *)fTrackHitsInfo->At(0);
 635   AliTrackHitsParam *param = (AliTrackHitsParam *)fTrackHitsParam->At(0);
 636   AliHitInfo * hinfo = (AliHitInfo *)fHitsPosAndQ->At(0,0);
 637
 638   if (!(info) || !(param) || !(hinfo) ) {
 639     fCurrentHit->fStatus = kFALSE;
 640     return kFALSE;
 641   }
 642
 643   fCurrentHit->fInfoIndex  = 0;
 644   fCurrentHit->fParamIndex = 0;
 645   fCurrentHit->fStackIndex = 0;
 646
 647   fCurrentHit->fHit.fSector = info->fVolumeID;
 648   fCurrentHit->fHit.SetTrack(info->fTrackID);
 649   fCurrentHit->fHit.SetX(param->fR*TMath::Cos(param->fFi));
 650   fCurrentHit->fHit.SetY(param->fR*TMath::Sin(param->fFi));
 651   fCurrentHit->fHit.SetZ(param->fZ);
 652   fCurrentHit->fHit.fQ = (Float_t)(hinfo->fCharge*AliTPCTrackHits::fgkTimePrecision);
 653   fCurrentHit->fHit.fTime = hinfo->fTime;
 654
 655   fCurrentHit->fR = param->fR;
 656
 657   return fCurrentHit->fStatus = kTRUE;
 658 }
 659
 660
 661 /*
 662 Bool_t AliTPCTrackHits::Next()
 663 {
 664   //
 665   //
 666   if (!(fCurrentHit->fStatus))
 667     return kFALSE;
 668
 669   fCurrentHit->fStackIndex++;
 670   AliHitInfo * hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex,
 671                                                       fCurrentHit->fStackIndex);
 672   AliTrackHitsInfo *info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
 673   AliTrackHitsParam *param =  (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
 674
 675   if (!hinfo) {
 676     hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex+1, 0);
 677     if (!hinfo)
 678       return fCurrentHit->fStatus = kFALSE;
 679     if (hinfo){
 680       fCurrentHit->fParamIndex++;
 681       fCurrentHit->fStackIndex = 0;
 682       param = (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
 683       if (!param)
 684         return fCurrentHit->fStatus = kFALSE;
 685       fCurrentHit->fR = param->fR;
 686
 687       if ((fCurrentHit->fInfoIndex+1<fTrackHitsInfo->GetSize())
 688         &&((info+1)->fHitParamIndex<=fCurrentHit->fParamIndex)){
 689         fCurrentHit->fInfoIndex++;
 690         info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
 691         if (!info)
 692           return fCurrentHit->fStatus = kFALSE;
 693         fCurrentHit->fHit.fSector = info->fVolumeID;
 694         fCurrentHit->fHit.SetTrack(info->fTrackID);
 695       }
 696     }
 697   }
 698   Double_t ratio;
 699   {
 700     //    Double_t dfi2 = param->fAn+2*param->fAd*(fCurrentHit->fR-param->fR);
 701     Double_t dfi2 = param->fAn;
 702     dfi2*=dfi2*fCurrentHit->fR*fCurrentHit->fR;
 703     //    Double_t ddz2 = param->fTheta+2*param->fThetaD*(fCurrentHit->fR-param->fR);
 704     Double_t ddz2 =  param->fTheta;
 705     ddz2*=ddz2;
 706     ratio = TMath::Sqrt(1.+ dfi2+ ddz2);
 707   }
 708
 709   fCurrentHit->fHit.fQ = hinfo->fCharge;
 710   fCurrentHit->fHit.fTime = (Float_t)(hinfo->fTime*AliTPCTrackHits::fgkTimePrecision);
 711   fCurrentHit->fR += fStep*hinfo->fHitDistance/ratio;
 712   Double_t dR = fCurrentHit->fR - param->fR;
 713   //Double_t dR =0;
 714   Double_t fi = param->fFi + (param->fAn*dR+param->fAd*dR*dR);
 715   Double_t z  = param->fZ + (param->fTheta*dR+param->fThetaD*dR*dR);
 716
 717   fCurrentHit->fHit.SetX(fCurrentHit->fR*TMath::Cos(fi));
 718   fCurrentHit->fHit.SetY(fCurrentHit->fR*TMath::Sin(fi));
 719   fCurrentHit->fHit.SetZ(z);
 720   return kTRUE;
 721 }
 722
 723 */
 724 AliTPChit * AliTPCTrackHits::GetHit() const
 725 {
 726   //
 727    return (fCurrentHit->fStatus)? &fCurrentHit->fHit:0;
 728   //return &fCurrentHit->fHit;
 729
 730 }
 731
 732
 733
 734 Bool_t AliTPCTrackHits::Next(Int_t id)
 735 {
 736   //
 737   //
 738   if (!(fCurrentHit->fStatus))
 739     return kFALSE;
 740
 741   //  fCurrentHit->fStackIndex++;
 742   AliHitInfo * hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex,
 743                                                       fCurrentHit->fStackIndex);
 744   AliTrackHitsInfo *info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
 745   if (!info) {
 746     fCurrentHit->fStatus = kFALSE;
 747     return kFALSE;
 748   }
 749   AliTrackHitsParam *param =  (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
 750
 751   if ( (id>=0) && (info!=0) && (info->fVolumeID!=id)){
 752     fCurrentHit->fInfoIndex++;
 753     info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
 754     if (!info) {
 755       fCurrentHit->fStatus = kFALSE;
 756       return kFALSE;
 757     }
 758     fCurrentHit->fParamIndex = info->fHitParamIndex;
 759     param =  (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
 760     fCurrentHit->fStackIndex =0;
 761     fCurrentHit->fR = param->fR;
 762     return Next(id);
 763   }
 764   if (!info) {
 765     fCurrentHit->fStatus = kFALSE;
 766     return kFALSE;
 767   }
 768   if (!hinfo) {
 769     hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex+1, 0);
 770     if (!hinfo){
 771       fCurrentHit->fStatus = kFALSE;
 772       return kFALSE;
 773     }
 774     if (hinfo){
 775       fCurrentHit->fParamIndex++;
 776       fCurrentHit->fStackIndex = 0;
 777       param = (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
 778       if (!param){
 779         fCurrentHit->fStatus = kFALSE;
 780         return kFALSE;
 781       }
 782       fCurrentHit->fR = param->fR;
 783
 784       if ((fCurrentHit->fInfoIndex+1<fTrackHitsInfo->GetSize())
 785         &&((info+1)->fHitParamIndex<=fCurrentHit->fParamIndex)){
 786         fCurrentHit->fInfoIndex++;
 787         info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
 788         if (!info){
 789           fCurrentHit->fStatus = kFALSE;
 790           return kFALSE;
 791         }
 792         if ( (id>=0) && (info!=0) && (info->fVolumeID!=id)){
 793           return Next(id);
 794         }
 795         fCurrentHit->fHit.fSector = info->fVolumeID;
 796         fCurrentHit->fHit.SetTrack(info->fTrackID);
 797       }
 798     }
 799   }
 800   Double_t ratio;
 801   {
 802     //    Double_t dfi2 = param->fAn+2*param->fAd*(fCurrentHit->fR-param->fR);
 803     Double_t dfi2 = param->fAn;
 804     dfi2*=dfi2*fCurrentHit->fR*fCurrentHit->fR;
 805     //    Double_t ddz2 = param->fTheta+2*param->fThetaD*(fCurrentHit->fR-param->fR);
 806     Double_t ddz2 =  param->fTheta;
 807     ddz2*=ddz2;
 808     ratio = TMath::Sqrt(1.+ dfi2+ ddz2);
 809   }
 810
 811   fCurrentHit->fHit.fQ = hinfo->fCharge;
 812   fCurrentHit->fHit.fTime = (Float_t)(hinfo->fTime*AliTPCTrackHits::fgkTimePrecision);
 813   fCurrentHit->fR += fStep*hinfo->fHitDistance/ratio;
 814   Double_t dR = fCurrentHit->fR - param->fR;
 815   //Double_t dR =0;
 816   Double_t fi = param->fFi + (param->fAn*dR+param->fAd*dR*dR);
 817   Double_t z  = param->fZ + (param->fTheta*dR+param->fThetaD*dR*dR);
 818
 819   fCurrentHit->fHit.SetX(fCurrentHit->fR*TMath::Cos(fi));
 820   fCurrentHit->fHit.SetY(fCurrentHit->fR*TMath::Sin(fi));
 821   fCurrentHit->fHit.SetZ(z);
 822   fCurrentHit->fHit.fSector = info->fVolumeID;
 823   fCurrentHit->fHit.SetTrack(info->fTrackID);
 824   //
 825   fCurrentHit->fStatus = kTRUE;
 826   fCurrentHit->fStackIndex++;
 827   return kTRUE;
 828 }
 829
 830
 831 AliTrackHitsParam * AliTPCTrackHits::GetParam()
 832 {
 833   //
 834   return (fCurrentHit->fStatus)? (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex) :0;
 835 }
 836
 837 AliHitInfo * AliTPCTrackHits::GetHitInfo()
 838 {
 839   //
 840   return (fCurrentHit->fStatus)?
 841     (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex,fCurrentHit->fStackIndex) :0;
 842 }
 843
 844
 845