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