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