Changes to obey the coding rules.
[u/mrichter/AliRoot.git] / TPC / AliTPCTrackHits.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 /*
17 $Log$
18 Revision 1.3  2000/11/02 10:22:50  kowal2
19 Logs added
20
21 */
22 ///////////////////////////////////////////////////////////////////////////////
23 //                                                                           //
24 //  Time Projection Chamber  track hits object                                //
25 //
26 //  Origin: Marian Ivanov , GSI Darmstadt
27 //
28 //  Class for storing simulated AliTPCHits  for given track                  //
29 //             -average compression comparing to classical ClonesArray is    //
30 //              around 5-7 (depending on the required hit precision)       //
31 //                                                                           //
32 //Begin_Html
33 /*
34 <img src="gif/AliTPCTrackHits.gif">
35 */
36 //End_Html
37 //                                                                           //
38 //                                                                          //
39 ///////////////////////////////////////////////////////////////////////////////
40
41 #include "TVector3.h"
42 #include "TClonesArray.h"    
43 #include "AliTPCTrackHits.h"
44 #include "AliTPC.h"
45
46 #include <iostream.h>
47
48
49
50 ClassImp(AliTPCTrackHits) 
51 LClassImp(AliTrackHitsInfo) 
52 LClassImp(AliTrackHitsParam)  
53 LClassImp(AliHitInfo)
54
55 Int_t AliTrackHitsInfo::fgCounter1 =0;
56 Int_t AliTrackHitsInfo::fgCounter2 =0;
57 Int_t AliTrackHitsParam::fgCounter1 =0;
58 Int_t AliTrackHitsParam::fgCounter2 =0;
59 Int_t AliHitInfo::fgCounter1 =0;
60 Int_t AliHitInfo::fgCounter2 =0;
61 Int_t AliTPCTrackHits::fgCounter1 =0;
62 Int_t AliTPCTrackHits::fgCounter2 =0;
63 const Double_t AliTPCTrackHits::fgkPrecision=1e-6;  //precision 
64 const Double_t AliTPCTrackHits::fgkPrecision2=1e-20;  //precision
65
66
67 /************************************************************/
68 //           Interface classes                              // 
69 #include "AliTPCTrackHitsInterfaces.h"
70
71
72
73
74 struct AliTPCCurrentHit {
75   AliTPChit fHit;
76   UInt_t   fInfoIndex;//   - current info pointer 
77   UInt_t   fParamIndex;//  - current param pointer
78   UInt_t   fStackIndex; // - current hit stack index
79   Double_t fR;   //current Radius
80   Bool_t  fStatus; //current status    
81 };   
82
83
84 struct  AliTPCTempHitInfo {
85   enum    { fkStackSize = 10000};
86   AliTPCTempHitInfo();   
87   void     NewParam(Double_t r, Double_t z, Double_t fi, Int_t q);
88   void     SetHit(Double_t r, Double_t z, Double_t fi, Int_t q);
89   Double_t * GetPosition(Int_t index){return &fPositionStack[index*3];}
90   void    UpdateParam(Double_t maxdelta); //recal
91   void   Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
92             Double_t fSumX,  Double_t fSumX2, Double_t fSumX3, 
93             Double_t fSumX4, Int_t n,
94               Double_t &a, Double_t &b, Double_t &c);
95   void  Fit(AliTrackHitsParam * param);
96   Double_t fSumDr;    //
97   Double_t fSumDr2;   //
98   Double_t fSumDr3;   //  
99   Double_t fSumDr4;   //
100   Double_t fSumDFi;  //
101   Double_t fSumDFiDr; //  
102   Double_t fSumDFiDr2;//
103   Double_t fSumDZ;     //
104   Double_t fSumDZDr;  //
105   Double_t fSumDZDr2;  //
106   Double_t fOldR;     //previos r
107   Double_t fPositionStack[3*fkStackSize];  //position stack 
108   UInt_t   fQStack[fkStackSize];           //Q stack
109   UInt_t fStackIndex;   //current stack index 
110   UInt_t fInfoIndex;    //current track info index
111   UInt_t fParamIndex;   //current track parameters index
112   AliTrackHitsInfo  * fInfo; //current track info
113   AliTrackHitsParam * fParam; //current track param
114 };
115
116
117 AliTPCTempHitInfo::AliTPCTempHitInfo()
118 {
119   //
120   //set to default value
121   fSumDr=fSumDr2=fSumDr3=fSumDr4=
122     fSumDFi=fSumDFiDr=fSumDFiDr2=
123     fSumDZ=fSumDZDr=fSumDZDr2=0;  
124   fStackIndex = 0;
125   fInfoIndex  = 0;
126   fParamIndex = 0;
127 }
128
129
130 void AliTPCTempHitInfo::NewParam(Double_t r, Double_t z, Double_t fi, Int_t q)
131 {
132   //
133   //reset stack and sum parameters
134   //store line initial point
135   fSumDr=fSumDr2=fSumDr3=fSumDr4=
136     fSumDFi=fSumDFiDr=fSumDFiDr2=
137     fSumDZ=fSumDZDr=fSumDZDr2=0;  
138   fStackIndex=0;
139   fParam->fR = r;
140   fOldR = r;
141   fParam->fZ = z;
142   fParam->fFi = fi;
143   fParam->fAn = 0.;
144   fParam->fAd = 0.;
145   fParam->fTheta =0.;
146   fParam->fThetaD =0.;
147   SetHit(r,z,fi,q);
148 }
149
150 void AliTPCTempHitInfo::SetHit(Double_t r, Double_t z, Double_t fi, Int_t q)
151 {
152   //
153   //add hit to the stack
154   //recalculate new estimete of line parameters
155   Double_t *f = GetPosition(fStackIndex);  
156   f[0] = r;
157   f[1] = z;
158   f[2] = fi;
159   fQStack[fStackIndex]=q;
160   if (fStackIndex==0) return;
161   Double_t dr  = (r-fParam->fR);
162   if (TMath::Abs(dr)<AliTPCTrackHits::fgkPrecision) dr =AliTPCTrackHits::fgkPrecision;
163   Double_t dfi = fi-fParam->fFi;
164   Double_t dz  = z -fParam->fZ; 
165   Double_t dr2 =dr*dr;
166   Double_t dr3 =dr2*dr;
167   Double_t dr4 =dr3*dr;
168   fSumDr +=dr;
169   fSumDr2+=dr2;
170   fSumDr3+=dr3;
171   fSumDr4+=dr4;
172   fSumDFi +=dfi;
173   fSumDFiDr+=dfi*dr;
174   fSumDFiDr2+=dfi*dr2;
175   fSumDZ +=dz;
176   fSumDZDr+=dz*dr;
177   fSumDZDr2+=dz*dr2;
178   
179   //update fit parameters
180   //
181   Double_t det = fSumDr2*fSumDr4-fSumDr3*fSumDr3;
182   if (TMath::Abs(det)<AliTPCTrackHits::fgkPrecision2) return;
183   if ( ( fStackIndex>1 )  ){
184     fParam->fAn = (fSumDr4*fSumDFiDr-fSumDr3*fSumDFiDr2)/det;
185     fParam->fAd = (fSumDr2*fSumDFiDr2-fSumDr3*fSumDFiDr)/det;
186   }
187   else
188     fParam->fAn = fSumDFiDr/fSumDr2;
189   if ( ( fStackIndex>1 )  ){
190     fParam->fTheta = (fSumDr4*fSumDZDr-fSumDr3*fSumDZDr2)/det;
191     fParam->fThetaD= (fSumDr2*fSumDZDr2-fSumDr3*fSumDZDr)/det;
192   }
193   else
194     fParam->fTheta = fSumDZDr/fSumDr2; 
195 }
196
197
198 void   AliTPCTempHitInfo::UpdateParam(Double_t maxdelta)
199 {
200   //recalc parameters not fixing origin point
201   if (fStackIndex>5){ 
202     Double_t a,b,c;
203     a=b=c=0;
204     Fit2(fSumDFi, fSumDFiDr, fSumDFiDr2, fSumDr,fSumDr2,fSumDr3,fSumDr4,
205          fStackIndex, a,b,c);
206     if (TMath::Abs(a)<maxdelta){
207       fParam->fFi +=a/fParam->fR;    
208       fParam->fAn = b;    
209       fParam->fAd = c;                  
210     }
211     Fit2(fSumDZ, fSumDZDr, fSumDZDr2, fSumDr,fSumDr2,fSumDr3,fSumDr4,
212          fStackIndex, a,b,c) ;   
213     if (TMath::Abs(a)<maxdelta){
214       fParam->fZ +=a;    
215       fParam->fTheta = b;    
216       fParam->fThetaD = c;   
217     }                         
218   }
219       
220 }
221 void   AliTPCTempHitInfo::Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
222             Double_t fSumX,  Double_t fSumX2, Double_t fSumX3, 
223             Double_t fSumX4, Int_t n,
224             Double_t &a, Double_t &b, Double_t &c)
225 {
226   //fit of second order
227   Double_t det = 
228     n* (fSumX2*fSumX4-fSumX3*fSumX3) -
229     fSumX*      (fSumX*fSumX4-fSumX3*fSumX2)+
230     fSumX2*     (fSumX*fSumX3-fSumX2*fSumX2);
231     
232   if (TMath::Abs(det)> AliTPCTrackHits::fgkPrecision) {    
233     a = 
234       (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
235        fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
236        fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det; 
237     b=
238       (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
239       fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
240       fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
241     c=
242       (n*(fSumX2*fSumYX2-fSumYX*fSumX3)-
243        fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
244        fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;  
245   }
246 }
247
248 void   AliTPCTempHitInfo::Fit(AliTrackHitsParam * param)
249 {
250   // fit fixing first and the last point 
251   //result stored in new param
252   Double_t dx2  = (GetPosition(fStackIndex))[0]-fParam->fR;
253   Double_t det = fSumDr4+dx2*fSumDr2-2*dx2*fSumDr3;
254   if ( (TMath::Abs(det)> AliTPCTrackHits::fgkPrecision) &&
255        ((TMath::Abs(dx2)> AliTPCTrackHits::fgkPrecision))){
256     Double_t dfi2 = (GetPosition(fStackIndex))[1]-fParam->fFi;
257     param->fAd = (fSumDFiDr2+dfi2*fSumDr-dx2*fSumDFiDr-dfi2*fSumDr3/dx2)/det;
258     param->fAn  = (dfi2-param->fAd*dx2*dx2)/dx2;
259     
260     Double_t dz2 = (GetPosition(fStackIndex))[1]-fParam->fZ;
261     param->fTheta = (fSumDZDr2+dz2*fSumDr-dx2*fSumDZDr-dz2*fSumDr3/dx2)/det;
262     param->fTheta  = (dz2-param->fAd*dx2*dx2)/dx2;
263   }
264   
265 }
266
267
268
269
270 AliTPCTrackHits::AliTPCTrackHits()
271 {
272   //
273   //default constructor
274   //
275   const Float_t kHitPrecision=0.002; //default precision for hit position in cm
276   const Float_t kStep =0.003;  //30 mum step 
277   const UShort_t kMaxDistance =100;  //maximum distance 100  
278
279   fPrecision=kHitPrecision; //precision in cm
280   fStep = kStep; //step size
281   fMaxDistance = kMaxDistance; //maximum distance
282   fTempInfo =0;
283   fTrackHitsInfo = new AliObjectArray("AliTrackHitsInfo"); 
284   fTrackHitsParam = new AliObjectArray("AliTrackHitsParam");
285   fHitsPosAndQ = new TArrayOfArrayVStack("AliHitInfo");
286   fCurrentHit = new AliTPCCurrentHit;
287   fgCounter1++;
288   fgCounter2++;
289
290
291
292 AliTPCTrackHits::~AliTPCTrackHits()
293 {
294   //
295   //default destructor
296   //
297   if (fTrackHitsInfo) delete fTrackHitsInfo;
298   if (fTrackHitsParam) delete fTrackHitsParam;
299   if (fHitsPosAndQ) delete fHitsPosAndQ;
300   if (fCurrentHit) delete fCurrentHit;
301   if (fTempInfo) delete fTempInfo;
302   fgCounter1--;
303 }
304
305 void AliTPCTrackHits::Clear()
306 {
307   //
308   //clear object 
309   fTrackHitsInfo->Clear();
310   fTrackHitsParam->Clear();  
311   //fTrackHitsInfo->Resize(0);
312   //fTrackHitsParam->Resize(0);
313   fHitsPosAndQ->Clear();
314
315   if (fTempInfo){
316     delete fTempInfo; 
317     fTempInfo =0;
318   } 
319 }
320
321
322 void AliTPCTrackHits::AddHitKartez(Int_t volumeID, Int_t trackID, Double_t x, 
323               Double_t y, Double_t z,Int_t q)
324 {
325   Double_t r = TMath::Sqrt(x*x+y*y);
326   Double_t fi = TMath::ACos(x/r);
327   if (y<0) fi*=-1.;
328     AddHit(volumeID,trackID,r,z,fi,q);
329 }
330
331 void AliTPCTrackHits::AddHit(Int_t volumeID, Int_t trackID, 
332                              Double_t r, Double_t z, Double_t fi, Int_t q)
333 {
334   //
335   Bool_t diff=kFALSE;
336   if (!fTempInfo) { //initialsation of track
337     fTempInfo = new AliTPCTempHitInfo;
338     //
339     if (fTrackHitsInfo->GetCapacity()<10) fTrackHitsInfo->Reserve(10);
340     fTrackHitsInfo->Resize(1);
341     fTempInfo->fInfoIndex =0;
342     if (fTrackHitsParam->GetCapacity()<100) fTrackHitsParam->Reserve(100);    
343     fTrackHitsParam->Resize(1);
344     //
345     fTempInfo->fInfo = 
346       (AliTrackHitsInfo*) (fTrackHitsInfo->At(0));
347     fTempInfo->fInfo->fVolumeID = volumeID;
348     fTempInfo->fInfo->fTrackID = trackID;
349     fTempInfo->fInfo->fHitParamIndex =0;     
350     fTempInfo->fInfoIndex = 0;
351     //
352     fTempInfo->fParam = 
353       (AliTrackHitsParam*) (fTrackHitsParam->At(0));
354     fTempInfo->fParamIndex = 0;
355     fTempInfo->NewParam(r,z,fi,q);
356     return;
357   }
358   
359   Int_t size = fHitsPosAndQ->GetSize();
360   if (size>(Int_t)fTempInfo->fParamIndex) {
361     fTempInfo->fParamIndex++;
362     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
363       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
364     fTempInfo->fParam = 
365       (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));  
366     fTempInfo->NewParam(r,z,fi,q);
367     return;
368   }
369   
370
371   // if new volume or new trackID  
372   if ( (volumeID!=fTempInfo->fInfo->fVolumeID) || 
373        (trackID!=fTempInfo->fInfo->fTrackID)){
374     diff=kTRUE;
375
376     FlushHitStack(kTRUE);
377       
378     fTempInfo->fInfoIndex++;
379     if (fTempInfo->fInfoIndex+1>fTrackHitsInfo->GetSize()) 
380       fTrackHitsInfo->Resize(fTempInfo->fInfoIndex+1);
381     fTempInfo->fInfo = 
382       (AliTrackHitsInfo*) (fTrackHitsInfo->At(fTempInfo->fInfoIndex));      
383     fTempInfo->fInfo->fVolumeID = volumeID;
384     fTempInfo->fInfo->fTrackID = trackID;
385     fTempInfo->fInfo->fHitParamIndex =fTempInfo->fParamIndex+1;  
386     // FlushHitStack(kTRUE);
387
388     fTempInfo->fParamIndex++;
389     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
390       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
391     fTempInfo->fParam = 
392       (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));  
393     fTempInfo->NewParam(r,z,fi,q);
394     return;
395   }
396      
397   //calculate current fit precission to next point
398   AliTrackHitsParam &param = *(fTempInfo->fParam);
399   Double_t dd=0;
400   Double_t dl=0;
401   Double_t ratio=0;
402   Double_t dr,dz,dfi,ddz,ddfi;
403   Double_t drhit,ddl;
404   dr=dz=dfi=ddz=ddfi=0;
405   drhit = r-fTempInfo->fOldR;
406   { 
407     //Double_t dfi2 = param.fAn+2*param.fAd*(r-param.fR); 
408     Double_t dfi2 = param.fAn;
409     dfi2*=dfi2*fTempInfo->fOldR*fTempInfo->fOldR;
410     //Double_t ddz2 =  param.fTheta+2*param.fThetaD*(r-param.fR);
411     Double_t ddz2 =  param.fTheta;
412     ddz2*=ddz2;
413     ratio = TMath::Sqrt(1.+ dfi2+ ddz2);  
414   }
415   dl = fStep * Short_t(TMath::Nint(drhit*ratio/fStep));
416   ddl = dl - drhit*ratio; 
417   fTempInfo->fOldR += dl/ratio; 
418
419   if (fTempInfo->fStackIndex>2){     
420     dr = r-param.fR;        
421     dz =  z-param.fZ;  
422     dfi = fi-param.fFi;
423     ddz = dr*param.fTheta+dr*dr*param.fThetaD-dz;
424     ddfi= dr*param.fAn+dr*dr*param.fAd-dfi;    
425     dd  = TMath::Sqrt(ddz*ddz+r*r*ddfi*ddfi+ddl*ddl); 
426     //
427   }        
428   //safety factor 1.25
429   if ( ( (dd*1.25>fPrecision) ) ||  
430        (fTempInfo->fStackIndex+4>fTempInfo->fkStackSize) || 
431        (TMath::Abs(dl/fStep)>fMaxDistance)  ) 
432     diff=kTRUE;
433   else{
434     fTempInfo->fStackIndex++;   
435     fTempInfo->SetHit(r,z,fi,q);
436     return;
437   }  
438   //if parameter changed 
439   if (FlushHitStack(kFALSE)){   //if full buffer flushed
440     fTempInfo->fParamIndex++;
441     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
442       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
443     fTempInfo->fParam = 
444       (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));  
445     fTempInfo->NewParam(r,z,fi,q);
446   }
447   else{
448     fTempInfo->fStackIndex++;
449     fTempInfo->SetHit(r,z,fi,q);              
450   }
451 }   
452
453 Bool_t AliTPCTrackHits::FlushHitStack(Bool_t force)
454 {
455   //
456   //write fHitsPosAndQ information from the stack to te arrays
457   if (!fTempInfo) return kFALSE; 
458   Int_t size = fHitsPosAndQ->GetSize();
459
460   if ( (size>0)&&(size!=(Int_t)fTempInfo->fParamIndex)) return kFALSE;  
461
462   if (fHitsPosAndQ->Push(fTempInfo->fStackIndex+1)!=fTempInfo->fParamIndex){
463     cout<<"internal error - contact MI\n";
464     return kFALSE;
465   }
466   AliHitInfo * info;
467  
468   AliTrackHitsParam & param = *(fTempInfo->fParam);
469   //recalculate track parameter not fixing first point
470   fTempInfo->UpdateParam(fStep/4.);
471   //fTempInfo->Fit(fTempInfo->fParam);  //- fixing the first and the last point
472
473   Double_t oldr = param.fR; 
474   //cout<<"C3"<<fTempInfo->fStackIndex<<"\n"<<flush;
475   UInt_t i;
476   Double_t dd;
477   for (i=0; i <= fTempInfo->fStackIndex; i++){
478     Double_t * position = fTempInfo->GetPosition(i);
479     Double_t   dr = position[0]-oldr;
480     Double_t   ratio; 
481     { 
482       //Double_t dfi2 = param.fAn+2*param.fAd*(position[0]-param.fR);
483       Double_t dfi2 = param.fAn;
484       dfi2*=dfi2*oldr*oldr;
485       //Double_t ddz2 =  param.fTheta+2*param.fThetaD*(position[0]-param.fR);
486       Double_t ddz2 =  param.fTheta;
487       ddz2*=ddz2;
488       ratio = TMath::Sqrt(1.+ dfi2+ ddz2);  
489     }
490
491     Double_t dl = fStep*(Short_t)TMath::Nint(dr*ratio/fStep);  
492     dr = dl/ratio; 
493     oldr+=dr;
494     //calculate precission
495     AliTrackHitsParam &param = *(fTempInfo->fParam);    
496     //real deltas
497     Double_t dr1=  position[0]-param.fR;
498     Double_t dz =  position[1]-param.fZ;
499
500     Double_t dfi = position[2]-param.fFi;
501     //extrapolated deltas
502     Double_t dr2 = oldr-param.fR; 
503     Double_t ddr = dr2-dr1;
504     Double_t ddz = dr2*param.fTheta+dr2*dr2*param.fThetaD-dz;
505     Double_t ddfi= dr2*param.fAn+dr2*dr2*param.fAd-dfi;    
506     dd = TMath::Sqrt(ddz*ddz+oldr*oldr*ddfi*ddfi+ddr*ddr); 
507
508     if ( (dd>fPrecision) ){ 
509       if (i==0){
510         param.fAn = 0;
511         param.fAd = 0;
512         param.fTheta =0;
513         param.fThetaD =0;
514         Double_t ddz = dr2*param.fTheta+dr2*dr2*param.fThetaD-dz;
515         Double_t ddfi= dr2*param.fAn+dr2*dr2*param.fAd-dfi;    
516         dl = 0;
517         dd = TMath::Sqrt(ddz*ddz+oldr*oldr*ddfi*ddfi+ddr*ddr); 
518       }
519       else
520         break;
521     }
522
523     info = (AliHitInfo*)(fHitsPosAndQ->At(fTempInfo->fParamIndex,i));
524     info->fHitDistance = Short_t(TMath::Nint(dl/fStep));
525     info->fCharge = Short_t(fTempInfo->fQStack[i]);
526     /*
527     cout<<"C2";
528     cout<<" "<<fTempInfo->fStackIndex<<" \t";
529     cout<<" "<<i<<" \t";
530     cout<<position[0]<<"\t";
531     cout<<position[1]<<"\t"; 
532     cout<<position[2]<<"\t";
533     cout<<param.fAn<<"\t";
534     cout<<param.fTheta<<"\t";
535     cout<<dr1<<"\t"<<ddr<<"\t"<<ddz<<"\t"<<ddfi<<"\t"<<dd<<"\n"<<flush;    
536     */
537   }    
538   
539   if (i<=fTempInfo->fStackIndex){ //if previous iteration not succesfull 
540     fHitsPosAndQ->Resize(fTempInfo->fParamIndex,i);
541     //
542     fTempInfo->fParamIndex++;
543     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
544       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
545     fTempInfo->fParam = 
546         (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));    
547     Double_t * p = fTempInfo->GetPosition(i);
548     UInt_t index2 = fTempInfo->fStackIndex;
549     fTempInfo->NewParam(p[0],p[1],p[2],fTempInfo->fQStack[i]);        
550     if (i+1<=index2) FlushHitStack2(i+1,index2);
551
552     if (force) return      FlushHitStack(kTRUE);      
553     return kFALSE;
554   }  
555   return kTRUE;
556
557  
558
559 void AliTPCTrackHits::FlushHitStack2(Int_t index1, Int_t index2)
560 {
561   //
562   // second iteration flush stack
563   // call only for hits where first iteration were not succesfully interpolated  
564   Double_t * positionstack = new Double_t[3*(index2-index1+1)];
565   UInt_t   * qstack        = new UInt_t[index2-index1+1];
566   memcpy(positionstack, &fTempInfo->fPositionStack[3*index1],
567          (3*(index2-index1+1))*sizeof(Double_t));
568   memcpy(qstack, &fTempInfo->fQStack[index1],(index2-index1+1)*sizeof(UInt_t));
569   Double_t *p = positionstack;
570   for (Int_t j=0; j<=index2-index1;j++){ 
571     fTempInfo->fStackIndex++;
572     fTempInfo->SetHit(p[3*j+0],p[3*j+1],p[3*j+2],qstack[j]);
573   }  
574   delete []positionstack;
575   delete []qstack;
576 }
577
578
579    
580
581
582
583
584   
585
586 Bool_t AliTPCTrackHits::First()
587 {
588   //
589   //set Current hit for the first hit
590   //
591   AliTrackHitsInfo *info = (AliTrackHitsInfo *)fTrackHitsInfo->At(0);
592   AliTrackHitsParam *param = (AliTrackHitsParam *)fTrackHitsParam->At(0);
593   AliHitInfo * hinfo = (AliHitInfo *)fHitsPosAndQ->At(0,0);
594
595   if (!(info) || !(param) || !(hinfo) ) {
596     fCurrentHit->fStatus = kFALSE;
597     return kFALSE;
598   }
599
600   fCurrentHit->fInfoIndex  = 0;
601   fCurrentHit->fParamIndex = 0;
602   fCurrentHit->fStackIndex = 0;
603
604   fCurrentHit->fHit.fSector = info->fVolumeID;
605   fCurrentHit->fHit.SetTrack(info->fTrackID);
606   fCurrentHit->fHit.SetX(param->fR*TMath::Cos(param->fFi));
607   fCurrentHit->fHit.SetY(param->fR*TMath::Sin(param->fFi));
608   fCurrentHit->fHit.SetZ(param->fZ); 
609   fCurrentHit->fHit.fQ = hinfo->fCharge;
610    
611   fCurrentHit->fR = param->fR;
612   
613   return fCurrentHit->fStatus = kTRUE;
614 }
615
616 Bool_t AliTPCTrackHits::Next()
617 {
618   //
619   //  
620   if (!(fCurrentHit->fStatus)) 
621     return kFALSE;
622
623   fCurrentHit->fStackIndex++;
624   AliHitInfo * hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex,
625                                                       fCurrentHit->fStackIndex);
626   AliTrackHitsInfo *info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
627   AliTrackHitsParam *param =  (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
628
629   if (!hinfo) {
630     hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex+1, 0);
631     if (!hinfo) 
632       return fCurrentHit->fStatus = kFALSE;
633     if (hinfo){ 
634       fCurrentHit->fParamIndex++;
635       fCurrentHit->fStackIndex = 0;
636       param = (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
637       if (!param) 
638         return fCurrentHit->fStatus = kFALSE;     
639       fCurrentHit->fR = param->fR;
640
641       if ((fCurrentHit->fInfoIndex+1<fTrackHitsInfo->GetSize())
642         &&((info+1)->fHitParamIndex<=fCurrentHit->fParamIndex)){
643         fCurrentHit->fInfoIndex++;
644         info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
645         if (!info) 
646           return fCurrentHit->fStatus = kFALSE;
647         fCurrentHit->fHit.fSector = info->fVolumeID;
648         fCurrentHit->fHit.SetTrack(info->fTrackID);
649       }
650     }  
651   } 
652   Double_t ratio;
653   { 
654     //    Double_t dfi2 = param->fAn+2*param->fAd*(fCurrentHit->fR-param->fR);
655     Double_t dfi2 = param->fAn;
656     dfi2*=dfi2*fCurrentHit->fR*fCurrentHit->fR;
657     //    Double_t ddz2 = param->fTheta+2*param->fThetaD*(fCurrentHit->fR-param->fR);
658     Double_t ddz2 =  param->fTheta;
659     ddz2*=ddz2;
660     ratio = TMath::Sqrt(1.+ dfi2+ ddz2);  
661   }
662
663   fCurrentHit->fHit.fQ = hinfo->fCharge;
664   fCurrentHit->fR += fStep*hinfo->fHitDistance/ratio;
665   Double_t dR = fCurrentHit->fR - param->fR;
666   //Double_t dR =0;
667   Double_t fi = param->fFi + (param->fAn*dR+param->fAd*dR*dR);
668   Double_t z  = param->fZ + (param->fTheta*dR+param->fThetaD*dR*dR);
669   
670   fCurrentHit->fHit.SetX(fCurrentHit->fR*TMath::Cos(fi));
671   fCurrentHit->fHit.SetY(fCurrentHit->fR*TMath::Sin(fi));
672   fCurrentHit->fHit.SetZ(z);  
673   return kTRUE;
674 }
675   
676 AliTPChit * AliTPCTrackHits::GetHit()
677 {
678   //
679    return (fCurrentHit->fStatus)? &fCurrentHit->fHit:0;
680   //return &fCurrentHit->fHit;
681
682 }  
683
684 AliTrackHitsParam * AliTPCTrackHits::GetParam()
685 {
686   //
687   return (fCurrentHit->fStatus)? (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex) :0;
688
689
690 AliHitInfo * AliTPCTrackHits::GetHitInfo()
691 {
692   //
693   return (fCurrentHit->fStatus)? 
694     (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex,fCurrentHit->fStackIndex) :0;
695
696
697
698