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