85666e987179728ebce41f6c2f762553d5c2bd08
[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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //  Time Projection Chamber  track hits object                                //
21 //
22 //  Origin: Marian Ivanov , GSI Darmstadt
23 //
24 //  Class for storing simulated AliTPCHits  for given track                  //
25 //             -average compression comparing to classical ClonesArray is    //
26 //              around 5-7 (depending on the required hit precision)       //
27 //                                                                           //
28 //Begin_Html
29 /*
30 <img src="gif/AliTPCTrackHits.gif">
31 */
32 //End_Html
33 //                                                                           //
34 //                                                                          //
35 ///////////////////////////////////////////////////////////////////////////////
36
37 //#include "TVector3.h"
38 //#include "TClonesArray.h"    
39 #include "AliTPCTrackHits.h"
40 #include "AliTPC.h"
41
42 //#include <Riostream.h>
43
44
45
46 ClassImp(AliTPCTrackHits) 
47 LClassImp(AliTrackHitsInfo) 
48 LClassImp(AliTrackHitsParam)  
49 LClassImp(AliHitInfo)
50
51 Int_t AliTrackHitsInfo::fgCounter1 =0;
52 Int_t AliTrackHitsInfo::fgCounter2 =0;
53 Int_t AliTrackHitsParam::fgCounter1 =0;
54 Int_t AliTrackHitsParam::fgCounter2 =0;
55 Int_t AliHitInfo::fgCounter1 =0;
56 Int_t AliHitInfo::fgCounter2 =0;
57 Int_t AliTPCTrackHits::fgCounter1 =0;
58 Int_t AliTPCTrackHits::fgCounter2 =0;
59 const Double_t AliTPCTrackHits::fgkPrecision=1e-6;  //precision 
60 const Double_t AliTPCTrackHits::fgkPrecision2=1e-20;  //precision
61
62
63 /************************************************************/
64 //           Interface classes                              // 
65 #include "AliTPCTrackHitsInterfaces.h"
66
67
68
69
70 struct AliTPCCurrentHit {
71   friend class AliTPCTrackHits;
72 private:
73   AliTPChit fHit;     //   - hit in "standard" representation
74   UInt_t   fInfoIndex;//   - current info pointer 
75   UInt_t   fParamIndex;//  - current param pointer
76   UInt_t   fStackIndex; // - current hit stack index
77   Double_t fR;   //current Radius
78   Bool_t  fStatus; //current status    
79 };   
80
81
82 struct  AliTPCTempHitInfo {
83   friend class AliTPCTrackHits;
84 private:
85   enum    { kStackSize = 100};
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;    //fSumDr
97   Double_t fSumDr2;   //fSumDr2
98   Double_t fSumDr3;   //  fSumDr3
99   Double_t fSumDr4;   //fSumDr4
100   Double_t fSumDFi;  //fSumDFi
101   Double_t fSumDFiDr; // fSumDFiDr 
102   Double_t fSumDFiDr2;//fSumDFiDr2
103   Double_t fSumDZ;     //fSumDZ
104   Double_t fSumDZDr;  //fSumDZDr
105   Double_t fSumDZDr2;  //fSumDZDr2
106   Double_t fOldR;     //previos r
107   Double_t fPositionStack[3*kStackSize];  //position stack 
108   UInt_t   fQStack[kStackSize];           //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(const AliTPCTrackHits& r) : TObject(r)
293 {
294 }
295
296
297 AliTPCTrackHits::~AliTPCTrackHits()
298 {
299   //
300   //default destructor
301   //
302   if (fTrackHitsInfo) delete fTrackHitsInfo;
303   if (fTrackHitsParam) delete fTrackHitsParam;
304   if (fHitsPosAndQ) delete fHitsPosAndQ;
305   if (fCurrentHit) delete fCurrentHit;
306   if (fTempInfo) delete fTempInfo;
307   fgCounter1--;
308 }
309
310 void AliTPCTrackHits::Clear()
311 {
312   //
313   //clear object 
314   fTrackHitsInfo->Clear();
315   fTrackHitsParam->Clear();  
316   //fTrackHitsInfo->Resize(0);
317   //fTrackHitsParam->Resize(0);
318   fHitsPosAndQ->Clear();
319
320   if (fTempInfo){
321     delete fTempInfo; 
322     fTempInfo =0;
323   } 
324 }
325
326
327 void AliTPCTrackHits::AddHitKartez(Int_t volumeID, Int_t trackID, Double_t x, 
328               Double_t y, Double_t z,Int_t q)
329 {
330   //add hits (cartesian)
331   Double_t r = TMath::Sqrt(x*x+y*y);
332   Double_t fi = TMath::ACos(x/r);
333   if (y<0) fi*=-1.;
334     AddHit(volumeID,trackID,r,z,fi,q);
335 }
336
337 void AliTPCTrackHits::AddHit(Int_t volumeID, Int_t trackID, 
338                              Double_t r, Double_t z, Double_t fi, Int_t q)
339 {
340   //
341   Bool_t diff=kFALSE;
342   if (!fTempInfo) { //initialsation of track
343     fTempInfo = new AliTPCTempHitInfo;
344     //
345     if (fTrackHitsInfo->GetCapacity()<10) fTrackHitsInfo->Reserve(10);
346     fTrackHitsInfo->Resize(1);
347     fTempInfo->fInfoIndex =0;
348     if (fTrackHitsParam->GetCapacity()<100) fTrackHitsParam->Reserve(100);    
349     fTrackHitsParam->Resize(1);
350     //
351     fTempInfo->fInfo = 
352       (AliTrackHitsInfo*) (fTrackHitsInfo->At(0));
353     fTempInfo->fInfo->fVolumeID = volumeID;
354     fTempInfo->fInfo->fTrackID = trackID;
355     fTempInfo->fInfo->fHitParamIndex =0;     
356     fTempInfo->fInfoIndex = 0;
357     //
358     fTempInfo->fParam = 
359       (AliTrackHitsParam*) (fTrackHitsParam->At(0));
360     fTempInfo->fParamIndex = 0;
361     fTempInfo->NewParam(r,z,fi,q);
362     return;
363   }
364   
365   Int_t size = fHitsPosAndQ->GetSize();
366   if (size>(Int_t)fTempInfo->fParamIndex) {
367     fTempInfo->fParamIndex++;
368     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
369       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
370     fTempInfo->fParam = 
371       (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));  
372     fTempInfo->NewParam(r,z,fi,q);
373     return;
374   }
375   
376
377   // if new volume or new trackID  
378   if ( (volumeID!=fTempInfo->fInfo->fVolumeID) || 
379        (trackID!=fTempInfo->fInfo->fTrackID)){
380     diff=kTRUE;
381
382     FlushHitStack(kTRUE);
383       
384     fTempInfo->fInfoIndex++;
385     if (fTempInfo->fInfoIndex+1>fTrackHitsInfo->GetSize()) 
386       fTrackHitsInfo->Resize(fTempInfo->fInfoIndex+1);
387     fTempInfo->fInfo = 
388       (AliTrackHitsInfo*) (fTrackHitsInfo->At(fTempInfo->fInfoIndex));      
389     fTempInfo->fInfo->fVolumeID = volumeID;
390     fTempInfo->fInfo->fTrackID = trackID;
391     fTempInfo->fInfo->fHitParamIndex =fTempInfo->fParamIndex+1;  
392     // FlushHitStack(kTRUE);
393
394     fTempInfo->fParamIndex++;
395     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
396       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
397     fTempInfo->fParam = 
398       (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));  
399     fTempInfo->NewParam(r,z,fi,q);
400     return;
401   }
402      
403   //calculate current fit precission to next point
404   AliTrackHitsParam &param = *(fTempInfo->fParam);
405   Double_t dd=0;
406   Double_t dl=0;
407   Double_t ratio=0;
408   Double_t dr,dz,dfi,ddz,ddfi;
409   Double_t drhit,ddl;
410   dr=dz=dfi=ddz=ddfi=0;
411   drhit = r-fTempInfo->fOldR;
412   { 
413     //Double_t dfi2 = param.fAn+2*param.fAd*(r-param.fR); 
414     Double_t dfi2 = param.fAn;
415     dfi2*=dfi2*fTempInfo->fOldR*fTempInfo->fOldR;
416     //Double_t ddz2 =  param.fTheta+2*param.fThetaD*(r-param.fR);
417     Double_t ddz2 =  param.fTheta;
418     ddz2*=ddz2;
419     ratio = TMath::Sqrt(1.+ dfi2+ ddz2);  
420   }
421   //
422  
423
424   dl = (TMath::Abs(drhit*ratio/fStep)<32000) ?  fStep * Short_t(TMath::Nint(drhit*ratio/fStep)):0;
425   ddl = dl - drhit*ratio; 
426   fTempInfo->fOldR += dl/ratio; 
427
428   if (fTempInfo->fStackIndex>2){     
429     dr = r-param.fR;        
430     dz =  z-param.fZ;  
431     dfi = fi-param.fFi;
432     ddz = dr*param.fTheta+dr*dr*param.fThetaD-dz;
433     ddfi= dr*param.fAn+dr*dr*param.fAd-dfi;    
434     dd  = TMath::Sqrt(ddz*ddz+r*r*ddfi*ddfi+ddl*ddl); 
435     //
436   }        
437   //safety factor 1.25
438   if ( ( (dd*1.25>fPrecision) ) ||  
439        (fTempInfo->fStackIndex+4>fTempInfo->kStackSize) || 
440        (TMath::Abs(dl/fStep)>fMaxDistance)  ) 
441     diff=kTRUE;
442   else{
443     fTempInfo->fStackIndex++;   
444     fTempInfo->SetHit(r,z,fi,q);
445     return;
446   }  
447   //if parameter changed 
448   if (FlushHitStack(kFALSE)){   //if full buffer flushed
449     fTempInfo->fParamIndex++;
450     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
451       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
452     fTempInfo->fParam = 
453       (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));  
454     fTempInfo->NewParam(r,z,fi,q);
455   }
456   else{
457     fTempInfo->fStackIndex++;
458     fTempInfo->SetHit(r,z,fi,q);              
459   }
460 }   
461
462 Bool_t AliTPCTrackHits::FlushHitStack(Bool_t force)
463 {
464   //
465   //write fHitsPosAndQ information from the stack to te arrays
466   if (!fTempInfo) return kFALSE; 
467   Int_t size = fHitsPosAndQ->GetSize();
468
469   if ( (size>0)&&(size!=(Int_t)fTempInfo->fParamIndex)) return kFALSE;  
470
471   if (fHitsPosAndQ->Push(fTempInfo->fStackIndex+1)!=fTempInfo->fParamIndex){
472     cout<<"internal error - contact MI\n";
473     return kFALSE;
474   }
475   AliHitInfo * info;
476  
477   AliTrackHitsParam & param = *(fTempInfo->fParam);
478   //recalculate track parameter not fixing first point
479   fTempInfo->UpdateParam(fStep/4.);
480   //fTempInfo->Fit(fTempInfo->fParam);  //- fixing the first and the last point
481
482   Double_t oldr = param.fR; 
483   //cout<<"C3"<<fTempInfo->fStackIndex<<"\n"<<flush;
484   UInt_t i;
485   Double_t dd;
486   for (i=0; i <= fTempInfo->fStackIndex; i++){
487     Double_t * position = fTempInfo->GetPosition(i);
488     Double_t   dr = position[0]-oldr;
489     Double_t   ratio; 
490     { 
491       //Double_t dfi2 = param.fAn+2*param.fAd*(position[0]-param.fR);
492       Double_t dfi2 = param.fAn;
493       dfi2*=dfi2*oldr*oldr;
494       //Double_t ddz2 =  param.fTheta+2*param.fThetaD*(position[0]-param.fR);
495       Double_t ddz2 =  param.fTheta;
496       ddz2*=ddz2;
497       ratio = TMath::Sqrt(1.+ dfi2+ ddz2);  
498     }
499
500     Double_t dl = (TMath::Abs(dr*ratio/fStep)<32000) ? fStep*(Short_t)TMath::Nint(dr*ratio/fStep):0;  
501     dr = dl/ratio; 
502     oldr+=dr;
503     //calculate precission
504     AliTrackHitsParam &param = *(fTempInfo->fParam);    
505     //real deltas
506     Double_t dr1=  position[0]-param.fR;
507     Double_t dz =  position[1]-param.fZ;
508
509     Double_t dfi = position[2]-param.fFi;
510     //extrapolated deltas
511     Double_t dr2 = oldr-param.fR; 
512     Double_t ddr = dr2-dr1;
513     Double_t ddz = dr2*param.fTheta+dr2*dr2*param.fThetaD-dz;
514     Double_t ddfi= dr2*param.fAn+dr2*dr2*param.fAd-dfi;    
515     dd = TMath::Sqrt(ddz*ddz+oldr*oldr*ddfi*ddfi+ddr*ddr); 
516
517     if ( (dd>fPrecision) ){ 
518       if (i==0){
519         param.fAn = 0;
520         param.fAd = 0;
521         param.fTheta =0;
522         param.fThetaD =0;
523         Double_t ddz = dr2*param.fTheta+dr2*dr2*param.fThetaD-dz;
524         Double_t ddfi= dr2*param.fAn+dr2*dr2*param.fAd-dfi;    
525         dl = 0;
526         dd = TMath::Sqrt(ddz*ddz+oldr*oldr*ddfi*ddfi+ddr*ddr); 
527       }
528       else
529         break;
530     }
531
532     info = (AliHitInfo*)(fHitsPosAndQ->At(fTempInfo->fParamIndex,i));
533     info->fHitDistance = (TMath::Abs(dl/fStep)<32000) ?Short_t(TMath::Nint(dl/fStep)):0;
534     info->fCharge = Short_t(fTempInfo->fQStack[i]);
535     /*
536     cout<<"C2";
537     cout<<" "<<fTempInfo->fStackIndex<<" \t";
538     cout<<" "<<i<<" \t";
539     cout<<position[0]<<"\t";
540     cout<<position[1]<<"\t"; 
541     cout<<position[2]<<"\t";
542     cout<<param.fAn<<"\t";
543     cout<<param.fTheta<<"\t";
544     cout<<dr1<<"\t"<<ddr<<"\t"<<ddz<<"\t"<<ddfi<<"\t"<<dd<<"\n"<<flush;    
545     */
546   }    
547   
548   if (i<=fTempInfo->fStackIndex){ //if previous iteration not succesfull 
549     fHitsPosAndQ->Resize(fTempInfo->fParamIndex,i);
550     //
551     fTempInfo->fParamIndex++;
552     if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
553       fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
554     fTempInfo->fParam = 
555         (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));    
556     Double_t * p = fTempInfo->GetPosition(i);
557     UInt_t index2 = fTempInfo->fStackIndex;
558     fTempInfo->NewParam(p[0],p[1],p[2],fTempInfo->fQStack[i]);        
559     if (i+1<=index2) FlushHitStack2(i+1,index2);
560
561     if (force) return      FlushHitStack(kTRUE);      
562     return kFALSE;
563   }  
564   return kTRUE;
565
566  
567
568 void AliTPCTrackHits::FlushHitStack2(Int_t index1, Int_t index2)
569 {
570   //
571   // second iteration flush stack
572   // call only for hits where first iteration were not succesfully interpolated  
573   Double_t * positionstack = new Double_t[3*(index2-index1+1)];
574   UInt_t   * qstack        = new UInt_t[index2-index1+1];
575   memcpy(positionstack, &fTempInfo->fPositionStack[3*index1],
576          (3*(index2-index1+1))*sizeof(Double_t));
577   memcpy(qstack, &fTempInfo->fQStack[index1],(index2-index1+1)*sizeof(UInt_t));
578   Double_t *p = positionstack;
579   for (Int_t j=0; j<=index2-index1;j++){ 
580     fTempInfo->fStackIndex++;
581     fTempInfo->SetHit(p[3*j+0],p[3*j+1],p[3*j+2],qstack[j]);
582   }  
583   delete []positionstack;
584   delete []qstack;
585 }
586
587
588    
589
590
591
592
593   
594
595 Bool_t AliTPCTrackHits::First()
596 {
597   //
598   //set Current hit for the first hit
599   //
600   AliTrackHitsInfo *info = (AliTrackHitsInfo *)fTrackHitsInfo->At(0);
601   AliTrackHitsParam *param = (AliTrackHitsParam *)fTrackHitsParam->At(0);
602   AliHitInfo * hinfo = (AliHitInfo *)fHitsPosAndQ->At(0,0);
603
604   if (!(info) || !(param) || !(hinfo) ) {
605     fCurrentHit->fStatus = kFALSE;
606     return kFALSE;
607   }
608
609   fCurrentHit->fInfoIndex  = 0;
610   fCurrentHit->fParamIndex = 0;
611   fCurrentHit->fStackIndex = 0;
612
613   fCurrentHit->fHit.fSector = info->fVolumeID;
614   fCurrentHit->fHit.SetTrack(info->fTrackID);
615   fCurrentHit->fHit.SetX(param->fR*TMath::Cos(param->fFi));
616   fCurrentHit->fHit.SetY(param->fR*TMath::Sin(param->fFi));
617   fCurrentHit->fHit.SetZ(param->fZ); 
618   fCurrentHit->fHit.fQ = hinfo->fCharge;
619    
620   fCurrentHit->fR = param->fR;
621   
622   return fCurrentHit->fStatus = kTRUE;
623 }
624
625
626 /*
627 Bool_t AliTPCTrackHits::Next()
628 {
629   //
630   //  
631   if (!(fCurrentHit->fStatus)) 
632     return kFALSE;
633
634   fCurrentHit->fStackIndex++;
635   AliHitInfo * hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex,
636                                                       fCurrentHit->fStackIndex);
637   AliTrackHitsInfo *info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
638   AliTrackHitsParam *param =  (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
639
640   if (!hinfo) {
641     hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex+1, 0);
642     if (!hinfo) 
643       return fCurrentHit->fStatus = kFALSE;
644     if (hinfo){ 
645       fCurrentHit->fParamIndex++;
646       fCurrentHit->fStackIndex = 0;
647       param = (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
648       if (!param) 
649         return fCurrentHit->fStatus = kFALSE;     
650       fCurrentHit->fR = param->fR;
651
652       if ((fCurrentHit->fInfoIndex+1<fTrackHitsInfo->GetSize())
653         &&((info+1)->fHitParamIndex<=fCurrentHit->fParamIndex)){
654         fCurrentHit->fInfoIndex++;
655         info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
656         if (!info) 
657           return fCurrentHit->fStatus = kFALSE;
658         fCurrentHit->fHit.fSector = info->fVolumeID;
659         fCurrentHit->fHit.SetTrack(info->fTrackID);
660       }
661     }  
662   } 
663   Double_t ratio;
664   { 
665     //    Double_t dfi2 = param->fAn+2*param->fAd*(fCurrentHit->fR-param->fR);
666     Double_t dfi2 = param->fAn;
667     dfi2*=dfi2*fCurrentHit->fR*fCurrentHit->fR;
668     //    Double_t ddz2 = param->fTheta+2*param->fThetaD*(fCurrentHit->fR-param->fR);
669     Double_t ddz2 =  param->fTheta;
670     ddz2*=ddz2;
671     ratio = TMath::Sqrt(1.+ dfi2+ ddz2);  
672   }
673
674   fCurrentHit->fHit.fQ = hinfo->fCharge;
675   fCurrentHit->fR += fStep*hinfo->fHitDistance/ratio;
676   Double_t dR = fCurrentHit->fR - param->fR;
677   //Double_t dR =0;
678   Double_t fi = param->fFi + (param->fAn*dR+param->fAd*dR*dR);
679   Double_t z  = param->fZ + (param->fTheta*dR+param->fThetaD*dR*dR);
680   
681   fCurrentHit->fHit.SetX(fCurrentHit->fR*TMath::Cos(fi));
682   fCurrentHit->fHit.SetY(fCurrentHit->fR*TMath::Sin(fi));
683   fCurrentHit->fHit.SetZ(z);  
684   return kTRUE;
685 }
686
687 */  
688 AliTPChit * AliTPCTrackHits::GetHit() const
689 {
690   //
691    return (fCurrentHit->fStatus)? &fCurrentHit->fHit:0;
692   //return &fCurrentHit->fHit;
693
694 }  
695
696
697
698 Bool_t AliTPCTrackHits::Next(Int_t id)
699 {
700   //
701   //  
702   if (!(fCurrentHit->fStatus)) 
703     return kFALSE;
704
705   //  fCurrentHit->fStackIndex++;
706   AliHitInfo * hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex,
707                                                       fCurrentHit->fStackIndex);
708   AliTrackHitsInfo *info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
709   if (!info) {
710     fCurrentHit->fStatus = kFALSE;
711     return kFALSE; 
712   }
713   AliTrackHitsParam *param =  (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
714
715   if ( (id>=0) && (info!=0) && (info->fVolumeID!=id)){
716     fCurrentHit->fInfoIndex++;
717     info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
718     if (!info) {
719       fCurrentHit->fStatus = kFALSE;
720       return kFALSE;
721     }
722     fCurrentHit->fParamIndex = info->fHitParamIndex;
723     param =  (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
724     fCurrentHit->fStackIndex =0;    
725     fCurrentHit->fR = param->fR;
726     return Next(id);
727   }
728   if (!info) {
729     fCurrentHit->fStatus = kFALSE;
730     return kFALSE; 
731   }
732   if (!hinfo) {
733     hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex+1, 0);
734     if (!hinfo){ 
735       fCurrentHit->fStatus = kFALSE;
736       return kFALSE;
737     }
738     if (hinfo){ 
739       fCurrentHit->fParamIndex++;
740       fCurrentHit->fStackIndex = 0;
741       param = (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
742       if (!param){ 
743         fCurrentHit->fStatus = kFALSE;
744         return kFALSE;     
745       }
746       fCurrentHit->fR = param->fR;
747
748       if ((fCurrentHit->fInfoIndex+1<fTrackHitsInfo->GetSize())
749         &&((info+1)->fHitParamIndex<=fCurrentHit->fParamIndex)){
750         fCurrentHit->fInfoIndex++;
751         info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
752         if (!info){ 
753           fCurrentHit->fStatus = kFALSE;
754           return kFALSE;
755         }
756         if ( (id>=0) && (info!=0) && (info->fVolumeID!=id)){
757           return Next(id);
758         }
759         fCurrentHit->fHit.fSector = info->fVolumeID;
760         fCurrentHit->fHit.SetTrack(info->fTrackID);
761       }
762     }  
763   } 
764   Double_t ratio;
765   { 
766     //    Double_t dfi2 = param->fAn+2*param->fAd*(fCurrentHit->fR-param->fR);
767     Double_t dfi2 = param->fAn;
768     dfi2*=dfi2*fCurrentHit->fR*fCurrentHit->fR;
769     //    Double_t ddz2 = param->fTheta+2*param->fThetaD*(fCurrentHit->fR-param->fR);
770     Double_t ddz2 =  param->fTheta;
771     ddz2*=ddz2;
772     ratio = TMath::Sqrt(1.+ dfi2+ ddz2);  
773   }
774
775   fCurrentHit->fHit.fQ = hinfo->fCharge;
776   fCurrentHit->fR += fStep*hinfo->fHitDistance/ratio;
777   Double_t dR = fCurrentHit->fR - param->fR;
778   //Double_t dR =0;
779   Double_t fi = param->fFi + (param->fAn*dR+param->fAd*dR*dR);
780   Double_t z  = param->fZ + (param->fTheta*dR+param->fThetaD*dR*dR);
781   
782   fCurrentHit->fHit.SetX(fCurrentHit->fR*TMath::Cos(fi));
783   fCurrentHit->fHit.SetY(fCurrentHit->fR*TMath::Sin(fi));
784   fCurrentHit->fHit.SetZ(z);  
785   fCurrentHit->fHit.fSector = info->fVolumeID;
786   fCurrentHit->fHit.SetTrack(info->fTrackID);
787   //
788   fCurrentHit->fStatus = kTRUE;
789   fCurrentHit->fStackIndex++;
790   return kTRUE;
791 }
792   
793
794 AliTrackHitsParam * AliTPCTrackHits::GetParam()
795 {
796   //
797   return (fCurrentHit->fStatus)? (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex) :0;
798
799
800 AliHitInfo * AliTPCTrackHits::GetHitInfo()
801 {
802   //
803   return (fCurrentHit->fStatus)? 
804     (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex,fCurrentHit->fStackIndex) :0;
805
806
807
808