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