Some of the coding violations corrected
[u/mrichter/AliRoot.git] / 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