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