Script supersceeded by AliForwarddNdetaTask.C and
[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 }
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 &paraml = *(fTempInfo->GetParam());    
602     //real deltas
603     Double_t dr1=  position[0]-paraml.GetR();
604     Double_t dz =  position[1]-paraml.GetZ();
605     Double_t dfi = position[2]-paraml.GetFi();
606     //extrapolated deltas
607     Double_t dr2 = oldr-paraml.GetR(); 
608     Double_t ddr = dr2-dr1;
609     Double_t ddz = dr2*paraml.GetTheta()+dr2*dr2*paraml.GetThetaD()-dz;
610     Double_t ddfi= dr2*paraml.GetAn()+dr2*dr2*paraml.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         paraml.SetAn(0);
618         paraml.SetAd(0);
619         paraml.SetTheta(0);
620         paraml.SetThetaD(0);
621         Double_t ddz1 = dr2*paraml.GetTheta()+dr2*dr2*paraml.GetThetaD()-dz;
622         Double_t ddfi1= dr2*paraml.GetAn()+dr2*dr2*paraml.GetAd()-dfi;    
623         dl = 0;
624         dd = TMath::Sqrt(ddz1*ddz1+oldr*oldr*ddfi1*ddfi1+ddr*ddr); 
625       }
626       else
627         break;
628     }
629
630     paraml.HitDistance(i)= Short_t(TMath::Nint(dl/fStep));
631     paraml.Charge(i)= Short_t(fTempInfo->GetQStack(i));
632     paraml.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