Small bug fixes to TPC files
[u/mrichter/AliRoot.git] / TPC / AliTPC.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.18  2000/04/17 09:37:33  kowal2
19 removed obsolete AliTPCDigitsDisplay.C
20
21 Revision 1.17.2.2  2000/04/10 08:15:12  kowal2
22
23 New, experimental data structure from M. Ivanov
24 New tracking algorithm
25 Different pad geometry for different sectors
26 Digitization rewritten
27
28 Revision 1.17.2.1  2000/04/10 07:56:53  kowal2
29 Not used anymore - removed
30
31 Revision 1.17  2000/01/19 17:17:30  fca
32 Introducing a list of lists of hits -- more hits allowed for detector now
33
34 Revision 1.16  1999/11/05 09:29:23  fca
35 Accept only signals > 0
36
37 Revision 1.15  1999/10/08 06:26:53  fca
38 Removed ClustersIndex - not used anymore
39
40 Revision 1.14  1999/09/29 09:24:33  fca
41 Introduction of the Copyright and cvs Log
42
43 */
44
45 ///////////////////////////////////////////////////////////////////////////////
46 //                                                                           //
47 //  Time Projection Chamber                                                  //
48 //  This class contains the basic functions for the Time Projection Chamber  //
49 //  detector. Functions specific to one particular geometry are              //
50 //  contained in the derived classes                                         //
51 //                                                                           //
52 //Begin_Html
53 /*
54 <img src="picts/AliTPCClass.gif">
55 */
56 //End_Html
57 //                                                                           //
58 //                                                                          //
59 ///////////////////////////////////////////////////////////////////////////////
60
61 #include <TMath.h>
62 #include <TRandom.h>
63 #include <TVector.h>
64 #include <TMatrix.h>
65 #include <TGeometry.h>
66 #include <TNode.h>
67 #include <TTUBS.h>
68 #include <TObjectTable.h>
69 #include "TParticle.h"
70 #include "AliTPC.h"
71 #include "AliRun.h"
72 #include <iostream.h>
73 #include <fstream.h>
74 #include "AliMC.h"
75
76
77 #include "AliTPCParam.h"
78 #include "AliTPCPRF2D.h"
79 #include "AliTPCRF1D.h"
80 #include "AliDigits.h"
81 #include "AliSimDigits.h"
82
83 #include "AliTPCDigitsArray.h"
84 #include "AliCluster.h"
85 #include "AliClusters.h"
86 #include "AliTPCClustersRow.h"
87 #include "AliTPCClustersArray.h"
88
89
90
91 ClassImp(AliTPC) 
92
93 //_____________________________________________________________________________
94 AliTPC::AliTPC()
95 {
96   //
97   // Default constructor
98   //
99   fIshunt   = 0;
100   fClusters = 0;
101   fHits     = 0;
102   fDigits   = 0;
103   fTracks   = 0;
104   fNsectors = 0;
105   fNtracks  = 0;
106   fNclusters= 0;
107   //MI changes
108   fDigitsArray = 0;
109   fClustersArray = 0;
110   fTPCParam = 0;
111 }
112  
113 //_____________________________________________________________________________
114 AliTPC::AliTPC(const char *name, const char *title)
115       : AliDetector(name,title)
116 {
117   //
118   // Standard constructor
119   //
120
121   //
122   // Initialise arrays of hits and digits 
123   fHits     = new TClonesArray("AliTPChit",  176);
124   gAlice->AddHitList(fHits);
125   //MI change  
126   fDigitsArray = 0;
127   fClustersArray= 0;
128   fTPCParam = 0;
129   //
130   // Initialise counters
131   fClusters = 0;
132   fTracks   = 0;
133   fNsectors = 0;
134   fNtracks  = 0;
135   fNclusters= 0;
136
137   //
138   fIshunt     =  0;
139   //
140   // Initialise color attributes
141   SetMarkerColor(kYellow);
142 }
143
144 //_____________________________________________________________________________
145 AliTPC::~AliTPC()
146 {
147   //
148   // TPC destructor
149   //
150   fIshunt   = 0;
151   delete fHits;
152   delete fDigits;
153   delete fClusters;
154   delete fTracks;
155   if (fDigitsArray!=0) delete fDigitsArray;
156   if (fClustersArray!=0) delete fClustersArray;
157
158   if (fTPCParam) delete fTPCParam;
159 }
160
161 //_____________________________________________________________________________
162 void AliTPC::AddCluster(Float_t *hits, Int_t *tracks)
163 {
164   //
165   // Add a simulated cluster to the list
166   //
167   if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000);
168   TClonesArray &lclusters = *fClusters;
169   new(lclusters[fNclusters++]) AliTPCcluster(hits,tracks);
170 }
171  
172 //_____________________________________________________________________________
173 void AliTPC::AddCluster(const AliTPCcluster &c)
174 {
175   //
176   // Add a simulated cluster copy to the list
177   //
178   if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",900000);
179   TClonesArray &lclusters = *fClusters;
180   new(lclusters[fNclusters++]) AliTPCcluster(c);
181 }
182  
183 //_____________________________________________________________________________
184 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
185 {
186   //
187   // Add a hit to the list
188   //
189   TClonesArray &lhits = *fHits;
190   new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
191 }
192  
193 //_____________________________________________________________________________
194 void AliTPC::AddTrack(Float_t *hits)
195 {
196   //
197   // Add a track to the list of tracks
198   //
199   TClonesArray &ltracks = *fTracks;
200   new(ltracks[fNtracks++]) AliTPCtrack(hits);
201 }
202
203 //_____________________________________________________________________________
204 void AliTPC::AddTrack(const AliTPCtrack& t)
205 {
206   //
207   // Add a track copy to the list of tracks
208   //
209   if(!fTracks) fTracks=new TClonesArray("AliTPCtrack",10000);
210   TClonesArray &ltracks = *fTracks;
211   new(ltracks[fNtracks++]) AliTPCtrack(t);
212 }
213
214 //_____________________________________________________________________________
215 void AliTPC::BuildGeometry()
216 {
217
218   //
219   // Build TPC ROOT TNode geometry for the event display
220   //
221   TNode *Node, *Top;
222   TTUBS *tubs;
223   Int_t i;
224   const int kColorTPC=19;
225   char name[5], title[25];
226   const Double_t kDegrad=TMath::Pi()/180;
227   const Double_t kRaddeg=180./TMath::Pi();
228
229
230   Float_t InnerOpenAngle = fTPCParam->GetInnerAngle();
231   Float_t OuterOpenAngle = fTPCParam->GetOuterAngle();
232
233   Float_t InnerAngleShift = fTPCParam->GetInnerAngleShift();
234   Float_t OuterAngleShift = fTPCParam->GetOuterAngleShift();
235
236   Int_t nLo = fTPCParam->GetNInnerSector()/2;
237   Int_t nHi = fTPCParam->GetNOuterSector()/2;  
238
239   const Double_t loAng = (Double_t)TMath::Nint(InnerOpenAngle*kRaddeg);
240   const Double_t hiAng = (Double_t)TMath::Nint(OuterOpenAngle*kRaddeg);
241   const Double_t loAngSh = (Double_t)TMath::Nint(InnerAngleShift*kRaddeg);
242   const Double_t hiAngSh = (Double_t)TMath::Nint(OuterAngleShift*kRaddeg);  
243
244
245   const Double_t loCorr = 1/TMath::Cos(0.5*loAng*kDegrad);
246   const Double_t hiCorr = 1/TMath::Cos(0.5*hiAng*kDegrad);
247
248   Double_t rl,ru;
249   
250
251   //
252   // Get ALICE top node
253   //
254
255   Top=gAlice->GetGeometry()->GetNode("alice");
256
257   //  inner sectors
258
259   rl = fTPCParam->GetInnerRadiusLow();
260   ru = fTPCParam->GetInnerRadiusUp();
261  
262
263   for(i=0;i<nLo;i++) {
264     sprintf(name,"LS%2.2d",i);
265     name[4]='\0';
266     sprintf(title,"TPC low sector %3d",i);
267     title[24]='\0';
268     
269     tubs = new TTUBS(name,title,"void",rl*loCorr,ru*loCorr,250.,
270                      loAng*(i-0.5)+loAngSh,loAng*(i+0.5)+loAngSh);
271     tubs->SetNumberOfDivisions(1);
272     Top->cd();
273     Node = new TNode(name,title,name,0,0,0,"");
274     Node->SetLineColor(kColorTPC);
275     fNodes->Add(Node);
276   }
277
278   // Outer sectors
279
280   rl = fTPCParam->GetOuterRadiusLow();
281   ru = fTPCParam->GetOuterRadiusUp();
282
283   for(i=0;i<nHi;i++) {
284     sprintf(name,"US%2.2d",i);
285     name[4]='\0';
286     sprintf(title,"TPC upper sector %d",i);
287     title[24]='\0';
288     tubs = new TTUBS(name,title,"void",rl*hiCorr,ru*hiCorr,250,
289                      hiAng*(i-0.5)+hiAngSh,hiAng*(i+0.5)+hiAngSh);
290     tubs->SetNumberOfDivisions(1);
291     Top->cd();
292     Node = new TNode(name,title,name,0,0,0,"");
293     Node->SetLineColor(kColorTPC);
294     fNodes->Add(Node);
295   }
296
297
298 }  
299   
300   
301
302 //_____________________________________________________________________________
303 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
304 {
305   //
306   // Calculate distance from TPC to mouse on the display
307   // Dummy procedure
308   //
309   return 9999;
310 }
311
312 //_____________________________________________________________________________
313 static Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
314 {
315   //
316   // Parametrised error of the cluster reconstruction (pad direction)   
317   //
318   pt=TMath::Abs(pt)*1000.;
319   Double_t x=r/pt;
320   tgl=TMath::Abs(tgl);
321   Double_t s=a_rphi - b_rphi*r*tgl + c_rphi*x*x + d_rphi*x;
322   if (s<0.4e-3) s=0.4e-3;
323   s*=1.3; //Iouri Belikov
324   return s;
325 }
326
327 //_____________________________________________________________________________
328 static Double_t SigmaZ2(Double_t r, Double_t tgl) 
329 {
330   //
331   // Parametrised error of the cluster reconstruction (drift direction)
332   //
333   tgl=TMath::Abs(tgl);
334   Double_t s=a_z - b_z*r*tgl + c_z*tgl*tgl;
335   if (s<0.4e-3) s=0.4e-3;
336   s*=1.3; //Iouri Belikov
337   return s;
338 }
339
340 //_____________________________________________________________________________
341 inline Double_t f1(Double_t x1,Double_t y1,
342                    Double_t x2,Double_t y2,
343                    Double_t x3,Double_t y3) 
344 {
345   //-----------------------------------------------------------------
346   // Initial approximation of the track curvature
347   //
348   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
349   //-----------------------------------------------------------------
350   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
351   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
352                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
353   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
354                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
355
356   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
357
358   return -xr*yr/sqrt(xr*xr+yr*yr); 
359 }
360
361
362 //_____________________________________________________________________________
363 inline Double_t f2(Double_t x1,Double_t y1,
364                    Double_t x2,Double_t y2,
365                    Double_t x3,Double_t y3) 
366 {
367   //-----------------------------------------------------------------
368   // Initial approximation of the track curvature times center of curvature
369   //
370   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
371   //-----------------------------------------------------------------
372   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
373   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
374                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
375   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
376                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
377
378   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
379   
380   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
381 }
382
383 //_____________________________________________________________________________
384 inline Double_t f3(Double_t x1,Double_t y1, 
385                    Double_t x2,Double_t y2,
386                    Double_t z1,Double_t z2) 
387 {
388   //-----------------------------------------------------------------
389   // Initial approximation of the tangent of the track dip angle
390   //
391   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
392   //-----------------------------------------------------------------
393   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
394 }
395
396 //_____________________________________________________________________________
397 static Int_t FindProlongation(AliTPCtrack& t, const AliTPCSector *sec,
398                             Int_t s, Int_t rf=0) 
399 {
400   //-----------------------------------------------------------------
401   // This function tries to find a track prolongation.
402   //
403   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
404   //-----------------------------------------------------------------
405   const Int_t ROWS_TO_SKIP=(t<10) ? 10 : Int_t(0.5*sec->GetNRows());
406   const Float_t MAX_CHI2=12.;
407   Int_t try_again=ROWS_TO_SKIP;
408   Double_t alpha=sec->GetAlpha();
409   Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5);
410
411   for (Int_t nr=sec->GetRowNumber(t.GetX())-1; nr>=rf; nr--) {
412     Double_t x=sec->GetX(nr), ymax=sec->GetMaxY(nr);
413     if (!t.PropagateTo(x)) return 0;
414
415     AliTPCcluster *cl=0;
416     Double_t max_chi2=MAX_CHI2;
417     const AliTPCRow& row=sec[s][nr];
418     Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
419     Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
420     Double_t road=5.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
421
422     if (road>30) {
423       if (t>4) cerr<<t<<" FindProlongation warning: Too broad road !\n"; 
424       return 0;
425     }
426
427     if (row) {
428       for (Int_t i=row.Find(y-road); i<row; i++) {
429         AliTPCcluster* c=(AliTPCcluster*)(row[i]);
430         if (c->fY > y+road) break;
431         if (c->IsUsed()) continue;
432         if ((c->fZ - z)*(c->fZ - z) > 25.*(t.GetSigmaZ2() + sz2)) continue;
433         Double_t chi2=t.GetPredictedChi2(c);
434         if (chi2 > max_chi2) continue;
435         max_chi2=chi2;
436         cl=c;       
437       }
438     }
439     if (cl) {
440       t.Update(cl,max_chi2);
441       cl->fdEdX=sec->GetPadPitchWidth()*TMath::Sqrt((1+t.GetTgl()*t.GetTgl())/
442                (1-(t.GetC()*x-t.GetEta())*(t.GetC()*x-t.GetEta())));
443       try_again=ROWS_TO_SKIP;
444     } else {
445       if (try_again==0) break;
446       if (y > ymax) {
447          s = (s+1) % ns;
448          if (!t.Rotate(alpha)) return 0;
449       } else if (y <-ymax) {
450          s = (s-1+ns) % ns;
451          if (!t.Rotate(-alpha)) return 0;
452       }
453       try_again--;
454     }
455   }
456
457   return 1;
458
459 }
460
461
462 //_____________________________________________________________________________
463 static void MakeSeeds(TObjArray& seeds,const AliTPCSector *sec, Int_t max_sec,
464 Int_t i1, Int_t i2)
465 {
466   //-----------------------------------------------------------------
467   // This function creates track seeds.
468   //
469   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
470   //-----------------------------------------------------------------
471   TMatrix C(5,5); TVector x(5);
472   Double_t alpha=sec->GetAlpha(), shift=sec->GetAlphaShift();
473   Double_t cs=cos(alpha), sn=sin(alpha);
474   for (Int_t ns=0; ns<max_sec; ns++) {
475     Int_t nl=sec[(ns-1+max_sec)%max_sec][i2];
476     Int_t nm=sec[ns][i2];
477     Int_t nu=sec[(ns+1)%max_sec][i2];
478     const AliTPCRow& r1=sec[ns][i1];
479     for (Int_t is=0; is < r1; is++) {
480       Double_t x1=sec->GetX(i1), y1=r1[is]->fY, z1=r1[is]->fZ;
481       for (Int_t js=0; js < nl+nm+nu; js++) {
482         const AliTPCcluster *cl;
483         Int_t ks;
484         Double_t x2=sec->GetX(i2), y2, z2, tmp;
485
486         if (js<nl) {
487           ks=(ns-1+max_sec)%max_sec;
488           const AliTPCRow& r2=sec[(ns-1+max_sec)%max_sec][i2];
489           cl=r2[js];
490           y2=cl->fY; z2=cl->fZ;
491           tmp= x2*cs+y2*sn;
492           y2 =-x2*sn+y2*cs; x2=tmp;
493         } else 
494           if (js<nl+nm) {
495             ks=ns;
496             const AliTPCRow& r2=sec[ns][i2];
497             cl=r2[js-nl];
498             y2=cl->fY; z2=cl->fZ;
499           } else {
500             ks=(ns+1)%max_sec;
501             const AliTPCRow& r2=sec[(ns+1)%max_sec][i2];
502             cl=r2[js-nl-nm];
503             y2=cl->fY; z2=cl->fZ;
504             tmp=x2*cs-y2*sn;
505             y2 =x2*sn+y2*cs; x2=tmp;
506           }
507
508         Double_t zz=z1 - z1/x1*(x1-x2); 
509         if (TMath::Abs(zz-z2)>5) continue;
510
511         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
512         if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
513
514         Double_t x3=0., y3=0.;//gRandom->Gaus(0.,TMath::Sqrt(cl->fSigmaY2));
515
516         x(0)=y1;
517         x(1)=z1;
518         x(2)=f1(x1,y1,x2,y2,x3,y3);
519         x(3)=f2(x1,y1,x2,y2,x3,y3);
520         x(4)=f3(x1,y1,x2,y2,z1,z2);
521         
522         if (TMath::Abs(x(2)*x1-x(3)) >= 0.999) continue;
523         
524         if (TMath::Abs(x(4)) > 1.2) continue;
525
526         Double_t a=asin(x(3));
527         Double_t zv=z1 - x(4)/x(2)*(a+asin(x(2)*x1-x(3)));
528         if (TMath::Abs(zv)>10.) continue; 
529
530         TMatrix X(6,6); X=0.; 
531         X(0,0)=r1[is]->fSigmaY2; X(1,1)=r1[is]->fSigmaZ2;
532         X(2,2)=cl->fSigmaY2;     X(3,3)=cl->fSigmaZ2;
533         X(4,4)=cl->fSigmaY2;     X(5,5)=cl->fSigmaZ2;
534         //X(4,4)=3./12.; X(5,5)=3./12.;
535         TMatrix F(5,6); F.UnitMatrix();
536         Double_t sy=sqrt(X(0,0)), sz=sqrt(X(1,1));
537         F(2,0)=(f1(x1,y1+sy,x2,y2,x3,y3)-x(2))/sy;
538         F(2,2)=(f1(x1,y1,x2,y2+sy,x3,y3)-x(2))/sy;
539         F(2,4)=(f1(x1,y1,x2,y2,x3,y3+sy)-x(2))/sy;
540         F(3,0)=(f2(x1,y1+sy,x2,y2,x3,y3)-x(3))/sy;
541         F(3,2)=(f2(x1,y1,x2,y2+sy,x3,y3)-x(3))/sy;
542         F(3,4)=(f2(x1,y1,x2,y2,x3,y3+sy)-x(3))/sy;
543         F(4,0)=(f3(x1,y1+sy,x2,y2,z1,z2)-x(4))/sy;
544         F(4,1)=(f3(x1,y1,x2,y2,z1+sz,z2)-x(4))/sz;
545         F(4,2)=(f3(x1,y1,x2,y2+sy,z1,z2)-x(4))/sy;
546         F(4,3)=(f3(x1,y1,x2,y2,z1,z2+sz)-x(4))/sz;
547         F(4,4)=0;
548         F(3,3)=0;
549         
550         TMatrix t(F,TMatrix::kMult,X);
551         C.Mult(t,TMatrix(TMatrix::kTransposed,F));
552
553         AliTPCtrack *track=new AliTPCtrack(r1[is], x, C, x1, ns*alpha+shift);
554         Int_t rc=FindProlongation(*track,sec,ns,i2);
555         if (rc<0 || *track<(i1-i2)/2) delete track;
556         else seeds.AddLast(track); 
557       }
558     }
559   }
560 }
561
562 //_____________________________________________________________________________
563 AliTPCParam *AliTPCSector::param;
564 void AliTPC::Clusters2Tracks()
565 {
566   //-----------------------------------------------------------------
567   // This is a track finder.
568   //
569   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
570   //-----------------------------------------------------------------
571   if (!fClusters) return;
572
573   AliTPCParam *p=fTPCParam;
574   AliTPCSector::SetParam(p);
575
576   const Int_t nis=p->GetNInnerSector()/2;
577   AliTPCSSector *ssec=new AliTPCSSector[nis];         
578   Int_t nrow_low=ssec->GetNRows();     
579
580   const Int_t nos=p->GetNOuterSector()/2;
581   AliTPCLSector *lsec=new AliTPCLSector[nos];
582   Int_t nrow_up=lsec->GetNRows();
583
584   Int_t ncl=fClusters->GetEntriesFast();
585   while (ncl--) {
586     AliTPCcluster *c=(AliTPCcluster*)fClusters->UncheckedAt(ncl);
587     Int_t sec=c->fSector, row=c->fPadRow;
588
589     if (sec<nis*2) {
590       ssec[sec%nis][row].InsertCluster(c);
591     } else {
592       sec -= nis*2;
593       lsec[sec%nos][row].InsertCluster(c);
594     }
595   }
596   
597   TObjArray seeds(20000);
598
599   Int_t nrows=nrow_low+nrow_up;
600   Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
601   MakeSeeds(seeds, lsec, nos, nrow_up-1, nrow_up-1-gap);
602   MakeSeeds(seeds, lsec, nos, nrow_up-1-shift, nrow_up-1-shift-gap);
603     
604   seeds.Sort();
605
606   Int_t found=0;
607   Int_t nseed=seeds.GetEntriesFast();
608
609   for (Int_t s=0; s<nseed; s++) {
610     AliTPCtrack *pt=(AliTPCtrack*)seeds.UncheckedAt(s), &t=*pt;
611     Double_t alpha=t.GetAlpha();
612     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
613     if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
614     Int_t ns=Int_t(alpha/lsec->GetAlpha())%nos;
615
616     if (!FindProlongation(t,lsec,ns)) continue;
617
618     alpha=t.GetAlpha() + 0.5*ssec->GetAlpha() - ssec->GetAlphaShift();
619     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
620     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
621     ns=Int_t(alpha/ssec->GetAlpha())%nis; //index of the inner sector needed
622
623     alpha=ns*ssec->GetAlpha() - t.GetAlpha();
624     if (!t.Rotate(alpha)) continue;
625     
626     if (!FindProlongation(t,ssec,ns)) continue;
627     
628     if (t >= Int_t(0.4*nrows)) {
629        AddTrack(t);
630        t.UseClusters();
631        cerr<<found++<<'\r';
632     }
633     delete pt; 
634   }  
635
636   delete[] ssec;
637   delete[] lsec;
638
639 }
640
641 //_____________________________________________________________________________
642 void AliTPC::CreateMaterials()
643 {
644   //-----------------------------------------------
645   // Create Materials for for TPC
646   //-----------------------------------------------
647
648   //-----------------------------------------------------------------
649   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
650   //-----------------------------------------------------------------
651
652   Int_t ISXFLD=gAlice->Field()->Integ();
653   Float_t SXMGMX=gAlice->Field()->Max();
654
655   Float_t amat[5]; // atomic numbers
656   Float_t zmat[5]; // z
657   Float_t wmat[5]; // proportions
658
659   Float_t density;
660
661   //  ********************* Gases *******************
662
663   //--------------------------------------------------------------
664   // pure gases
665   //--------------------------------------------------------------
666
667   // Ne
668
669
670   Float_t a_ne = 20.18;
671   Float_t z_ne = 10.;
672   
673   density = 0.0009;
674
675   AliMaterial(20,"Ne",a_ne,z_ne,density,999.,999.);
676
677   // Ar
678
679   Float_t a_ar = 39.948;
680   Float_t z_ar = 18.;
681
682   density = 0.001782;
683  
684   AliMaterial(21,"Ar",a_ar,z_ar,density,999.,999.);
685
686   Float_t a_pure[2];
687   
688   a_pure[0] = a_ne;
689   a_pure[1] = a_ar;
690   
691
692   //--------------------------------------------------------------
693   // gases - compounds
694   //--------------------------------------------------------------
695
696   Float_t amol[3];
697
698   //  CO2
699
700   amat[0]=12.011;
701   amat[1]=15.9994;
702
703   zmat[0]=6.;
704   zmat[1]=8.;
705
706   wmat[0]=1.;
707   wmat[1]=2.;
708
709   density=0.001977;
710
711   amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
712
713   AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
714
715   // CF4
716
717   amat[0]=12.011;
718   amat[1]=18.998;
719
720   zmat[0]=6.;
721   zmat[1]=9.;
722  
723   wmat[0]=1.;
724   wmat[1]=4.;
725  
726   density=0.003034;
727
728   amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
729
730   AliMixture(11,"CF4",amat,zmat,density,-2,wmat); 
731
732   // CH4
733
734   amat[0]=12.011;
735   amat[1]=1.;
736
737   zmat[0]=6.;
738   zmat[1]=1.;
739
740   wmat[0]=1.;
741   wmat[1]=4.;
742
743   density=0.000717;
744
745   amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
746
747   AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
748
749   //----------------------------------------------------------------
750   // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
751   //----------------------------------------------------------------
752  
753   char namate[21];
754  
755   density = 0.;
756   Float_t am=0;
757   Int_t nc;
758
759   Float_t a,z,rho,absl,X0,buf[1];
760   Int_t nbuf;
761
762   for(nc = 0;nc<fNoComp;nc++)
763     {
764     
765       // retrive material constants
766       
767       gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf);
768
769       amat[nc] = a;
770       zmat[nc] = z;
771
772       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
773  
774       am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? a_pure[nnc] : amol[nnc]); 
775       density += fMixtProp[nc]*rho;  // density of the mixture
776       
777     }
778
779   // mixture proportions by weight!
780
781   for(nc = 0;nc<fNoComp;nc++)
782     {
783
784       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
785
786       wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ? a_pure[nnc] : amol[nnc])/am;
787
788     }  
789   
790   AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
791   AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
792   AliMixture(33,"Drift gas 3",amat,zmat,density,fNoComp,wmat); 
793
794   AliMedium(2, "Drift gas 1", 31, 0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
795   AliMedium(3, "Drift gas 2", 32, 0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
796   AliMedium(4, "Drift gas 3", 33, 1, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
797
798   // Air 
799
800   AliMaterial(24, "Air", 14.61, 7.3, .001205, 30420., 67500.);
801
802   AliMedium(24, "Air", 24, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
803
804   //----------------------------------------------------------------------
805   //               solid materials
806   //----------------------------------------------------------------------
807
808   // Al
809
810   AliMaterial(30, "Al", 26.98, 13., 2.7, 8.9, 37.2);
811
812   AliMedium(0, "Al",30, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1);
813
814   // Si
815
816   AliMaterial(31, "Si", 28.086, 14.,2.33, 9.36, 999.);
817
818   AliMedium(7, "Al",31, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1);
819   
820
821   // Mylar C5H4O2
822
823   amat[0]=12.011;
824   amat[1]=1.;
825   amat[2]=15.9994;
826
827   zmat[0]=6.;
828   zmat[1]=1.;
829   zmat[2]=8.;
830
831   wmat[0]=5.;
832   wmat[1]=4.;
833   wmat[2]=2.; 
834
835   density = 1.39;
836   
837   AliMixture(32, "Mylar",amat,zmat,density,-3,wmat);
838
839   AliMedium(5, "Mylar",32, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
840
841
842
843
844   // Carbon (normal)
845
846   AliMaterial(33,"C normal",12.011,6.,2.265,18.8,999.);
847
848   AliMedium(6,"C normal",33,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
849
850   // G10 for inner and outr field cage
851   // G10 is 60% SiO2 + 40% epoxy, right now I use A and Z for SiO2
852
853   Float_t rhoFactor;
854
855   amat[0]=28.086;
856   amat[1]=15.9994;
857
858   zmat[0]=14.;
859   zmat[1]=8.;
860
861   wmat[0]=1.;
862   wmat[1]=2.;
863
864   density = 1.7;
865   
866
867   AliMixture(34,"G10 aux.",amat,zmat,density,-2,wmat);
868
869
870   gMC->Gfmate((*fIdmate)[34],namate,a,z,rho,X0,absl,buf,nbuf);
871
872   Float_t thickX0 = 0.0052; // field cage in X0 units
873   
874   Float_t thick = 2.; // in cm
875
876   X0=19.4; // G10 
877
878   rhoFactor = X0*thickX0/thick;
879   density = rho*rhoFactor;
880
881   AliMaterial(35,"G10-fc",a,z,density,999.,999.);
882
883   AliMedium(8,"G10-fc",35,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
884
885   thickX0 = 0.0027; // inner vessel (eta <0.9)
886   thick=0.5;
887   rhoFactor = X0*thickX0/thick;
888   density = rho*rhoFactor;
889
890   AliMaterial(36,"G10-iv",a,z,density,999.,999.);  
891
892   AliMedium(9,"G10-iv",36,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
893
894   //  Carbon fibre  
895   
896   gMC->Gfmate((*fIdmate)[33],namate,a,z,rho,X0,absl,buf,nbuf);
897
898   thickX0 = 0.0133; // outer vessel
899   thick=3.0;
900   rhoFactor = X0*thickX0/thick;
901   density = rho*rhoFactor;
902
903
904   AliMaterial(37,"C-ov",a,z,density,999.,999.);
905
906   AliMedium(10,"C-ov",37,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);  
907
908   thickX0=0.015; // inner vessel (cone, eta > 0.9)
909   thick=1.5;
910   rhoFactor = X0*thickX0/thick;
911   density = rho*rhoFactor;
912
913   AliMaterial(38,"C-ivc",a,z,density,999.,999.);
914
915   AliMedium(11,"C-ivc",38,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01);
916
917   //
918
919   AliMedium(12,"CO2",10,0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001);
920     
921 }
922
923 //_____________________________________________________________________________
924 struct Bin {
925    UShort_t q;
926    UInt_t mask;
927    Bin();
928 };
929 Bin::Bin() {q=0; mask=0xFFFFFFFE;}
930
931 struct Peak {
932    Int_t k;
933    UInt_t mask;
934 };
935 inline Bool_t IsMaximum(Int_t k, Int_t max, const Bin *bins) {
936   UShort_t q=bins[k].q;
937   if (q==1023) return kFALSE;
938   if (bins[k-max].q > q) return kFALSE;
939   if (bins[k-1  ].q > q) return kFALSE; 
940   if (bins[k+max].q > q) return kFALSE; 
941   if (bins[k+1  ].q > q) return kFALSE; 
942   if (bins[k-max-1].q > q) return kFALSE;
943   if (bins[k+max-1].q > q) return kFALSE; 
944   if (bins[k+max+1].q > q) return kFALSE; 
945   if (bins[k-max+1].q > q) return kFALSE;
946   return kTRUE; 
947 }
948 static void FindPeaks(Int_t k, Int_t max, Bin *bins, Peak *peaks, Int_t& n) {
949 //if (n>=31) return;
950   if (n<31)
951   if (IsMaximum(k,max,bins)) {
952     peaks[n].k=k; peaks[n].mask=(2<<n);
953     n++;
954   }
955   bins[k].mask=0;
956   if (bins[k-max].mask&1) FindPeaks(k-max,max,bins,peaks,n);
957   if (bins[k-1  ].mask&1) FindPeaks(k-1  ,max,bins,peaks,n);
958   if (bins[k+max].mask&1) FindPeaks(k+max,max,bins,peaks,n);
959   if (bins[k+1  ].mask&1) FindPeaks(k+1  ,max,bins,peaks,n);
960 }
961
962 static void MarkPeak(Int_t k, Int_t max, Bin *bins, UInt_t m) {
963   UShort_t q=bins[k].q;
964
965   bins[k].mask |= m; 
966
967   if (bins[k-max].q <= q)
968      if ((bins[k-max].mask&m) == 0) MarkPeak(k-max,max,bins,m);
969   if (bins[k-1  ].q <= q)
970      if ((bins[k-1  ].mask&m) == 0) MarkPeak(k-1  ,max,bins,m);
971   if (bins[k+max].q <= q)
972      if ((bins[k+max].mask&m) == 0) MarkPeak(k+max,max,bins,m);
973   if (bins[k+1  ].q <= q)
974      if ((bins[k+1  ].mask&m) == 0) MarkPeak(k+1  ,max,bins,m);
975 }
976
977 static void MakeCluster(Int_t k,Int_t max,Bin *bins,UInt_t m,AliTPCcluster &c){
978   Float_t q=(Float_t)bins[k].q;
979   Int_t i=k/max, j=k-i*max;
980   c.fY += i*q;
981   c.fZ += j*q;
982   c.fSigmaY2 += i*i*q;
983   c.fSigmaZ2 += j*j*q;
984   c.fQ += q;
985
986   bins[k].mask = 0xFFFFFFFE;
987   
988   if (bins[k-max].mask == m) MakeCluster(k-max,max,bins,m,c);
989   if (bins[k-1  ].mask == m) MakeCluster(k-1  ,max,bins,m,c);
990   if (bins[k+max].mask == m) MakeCluster(k+max,max,bins,m,c);
991   if (bins[k+1  ].mask == m) MakeCluster(k+1  ,max,bins,m,c);
992 }
993
994 //_____________________________________________________________________________
995 void AliTPC::Digits2Clusters()
996 {
997   //-----------------------------------------------------------------
998   // This is a simple cluster finder.
999   //
1000   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
1001   //-----------------------------------------------------------------
1002   AliTPCParam *par = fTPCParam;
1003   const Int_t MAXZ=par->GetMaxTBin()+2;
1004
1005   TTree *t = (TTree *)gDirectory->Get("TreeD_75x40_100x60");
1006   AliSimDigits digarr, *dummy=&digarr;
1007   t->GetBranch("Segment")->SetAddress(&dummy);
1008   Stat_t sectors_by_rows = t->GetEntries();
1009   for (Int_t n=0; n<sectors_by_rows; n++) {
1010     t->GetEvent(n);
1011     Int_t sec, row;
1012     if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
1013        cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
1014        continue;
1015     }
1016
1017     Float_t rx=par->GetPadRowRadii(sec,row);
1018
1019     Int_t npads, sign;
1020     {
1021        Int_t nis=par->GetNInnerSector(), nos=par->GetNOuterSector();
1022        if (sec < nis) {
1023           npads = par->GetNPadsLow(row);
1024           sign = (sec < nis/2) ? 1 : -1;
1025        } else {
1026           npads = par->GetNPadsUp(row);
1027           sign = ((sec-nis) < nos/2) ? 1 : -1;
1028        }
1029     }
1030
1031     const Int_t MAXBIN=MAXZ*(npads+2);
1032     Bin *bins=new Bin[MAXBIN];
1033
1034     digarr.First();
1035     do {
1036        Short_t dig=digarr.CurrentDigit();
1037        if (dig<=par->GetZeroSup()) continue;
1038        Int_t j=digarr.CurrentRow()+1, i=digarr.CurrentColumn()+1;
1039        bins[i*MAXZ+j].q=dig;
1040        bins[i*MAXZ+j].mask=1;
1041     } while (digarr.Next());
1042
1043     Int_t ncl=0;
1044     for (Int_t i=0; i<MAXBIN; i++) {
1045       if ((bins[i].mask&1) == 0) continue;
1046       Peak peaks[32]; Int_t npeaks=0;
1047       FindPeaks(i, MAXZ, bins, peaks, npeaks);
1048
1049       if (npeaks>30) continue;
1050
1051       Int_t k,l;
1052       for (k=0; k<npeaks-1; k++){//mark adjacent peaks
1053         if (peaks[k].k < 0) continue; //this peak is already removed
1054         for (l=k+1; l<npeaks; l++) {
1055            if (peaks[l].k < 0) continue; //this peak is already removed
1056            Int_t ki=peaks[k].k/MAXZ, kj=peaks[k].k - ki*MAXZ;
1057            Int_t li=peaks[l].k/MAXZ, lj=peaks[l].k - li*MAXZ;
1058            Int_t di=TMath::Abs(ki - li);
1059            Int_t dj=TMath::Abs(kj - lj);
1060            if (di>1 || dj>1) continue;
1061            if (bins[peaks[k].k].q > bins[peaks[l].k].q) {
1062               peaks[l].mask=peaks[k].mask;
1063               peaks[l].k*=-1;
1064            } else {
1065               peaks[k].mask=peaks[l].mask;
1066               peaks[k].k*=-1;
1067               break;
1068            } 
1069         }
1070       }
1071
1072       for (k=0; k<npeaks; k++) {
1073         MarkPeak(TMath::Abs(peaks[k].k), MAXZ, bins, peaks[k].mask);
1074       }
1075         
1076       for (k=0; k<npeaks; k++) {
1077          if (peaks[k].k < 0) continue; //removed peak
1078          AliTPCcluster c;
1079          MakeCluster(peaks[k].k, MAXZ, bins, peaks[k].mask, c);
1080          if (c.fQ < 5) continue; //noise cluster
1081          c.fY /= c.fQ;
1082          c.fZ /= c.fQ;
1083
1084          Double_t s2 = c.fSigmaY2/c.fQ - c.fY*c.fY;
1085          c.fSigmaY2 = s2 + 1./12.;
1086          c.fSigmaY2 *= par->GetPadPitchWidth(sec)*par->GetPadPitchWidth(sec);
1087          if (s2 != 0.) {
1088             c.fSigmaY2 *= 0.064*1.3*1.3;
1089             if (sec<par->GetNInnerSector()) c.fSigmaY2 *= 1.44*1.44;
1090          }
1091
1092          s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ;
1093          c.fSigmaZ2 = s2 + 1./12.;
1094          c.fSigmaZ2 *= par->GetZWidth()*par->GetZWidth();
1095          if (s2 != 0.) {
1096             c.fSigmaZ2 *= 0.10*1.3*1.3;
1097             if (sec<par->GetNInnerSector()) c.fSigmaZ2 *= 1.33*1.33;
1098          }
1099
1100          c.fY = (c.fY - 0.5 - 0.5*npads)*par->GetPadPitchWidth(sec);
1101          c.fZ = par->GetZWidth()*(c.fZ-1); 
1102          c.fZ -= 3.*par->GetZSigma(); // PASA delay 
1103          c.fZ = sign*(z_end - c.fZ);
1104
1105          if (rx<230./250.*TMath::Abs(c.fZ)) continue;
1106
1107          c.fSector=sec;
1108          c.fPadRow=row;
1109          Int_t ki=peaks[k].k/MAXZ, kj=peaks[k].k - ki*MAXZ;
1110          c.fTracks[0]=digarr.GetTrackID(kj-1,ki-1,0);
1111          c.fTracks[1]=digarr.GetTrackID(kj-1,ki-1,1);
1112          c.fTracks[2]=digarr.GetTrackID(kj-1,ki-1,2);
1113
1114          c.fQ=bins[peaks[k].k].q;
1115
1116          if (ki==1 || ki==npads || kj==1 || kj==MAXZ-2) {
1117            c.fSigmaY2 *= 25.;
1118            c.fSigmaZ2 *= 4.;
1119          }
1120
1121          AddCluster(c); ncl++;
1122       }
1123     }
1124
1125     cerr<<"sector, row, compressed digits, clusters: "
1126     <<sec<<' '<<row<<' '<<digarr.GetSize()<<' '<<ncl<<"                  \r";
1127         
1128     delete[] bins;  
1129   }
1130 }
1131
1132 //_____________________________________________________________________________
1133 void AliTPC::Hits2Clusters()
1134 {
1135   //--------------------------------------------------------
1136   // TPC simple cluster generator from hits
1137   // obtained from the TPC Fast Simulator
1138   // The point errors are taken from the parametrization
1139   //--------------------------------------------------------
1140
1141   //-----------------------------------------------------------------
1142   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1143   //-----------------------------------------------------------------
1144
1145    if(fTPCParam == 0){
1146      printf("AliTPCParam MUST be created firstly\n");
1147      return;
1148    }
1149
1150   Float_t sigma_rphi,sigma_z,cl_rphi,cl_z;
1151   //
1152   TParticle *particle; // pointer to a given particle
1153   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1154   TClonesArray *Particles; //pointer to the particle list
1155   Int_t sector,nhits;
1156   Int_t ipart;
1157   Float_t xyz[5];
1158   Float_t pl,pt,tanth,rpad,ratio;
1159   Float_t cph,sph;
1160   
1161   //---------------------------------------------------------------
1162   //  Get the access to the tracks 
1163   //---------------------------------------------------------------
1164   
1165   TTree *TH = gAlice->TreeH();
1166   Stat_t ntracks = TH->GetEntries();
1167   
1168   //------------------------------------------------------------
1169   // Loop over all sectors (72 sectors for 20 deg
1170   // segmentation for both lower and upper sectors)
1171   // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
1172   // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
1173   //
1174   // First cluster for sector 0 starts at "0"
1175   //------------------------------------------------------------
1176    
1177   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++){
1178     //MI change
1179     fTPCParam->AdjustCosSin(isec,cph,sph);
1180     
1181     //------------------------------------------------------------
1182     // Loop over tracks
1183     //------------------------------------------------------------
1184     
1185     for(Int_t track=0;track<ntracks;track++){
1186       ResetHits();
1187       TH->GetEvent(track);
1188       //
1189       //  Get number of the TPC hits and a pointer
1190       //  to the particles
1191       //
1192       nhits=fHits->GetEntriesFast();
1193       Particles=gAlice->Particles();
1194       //
1195       // Loop over hits
1196       //
1197       for(Int_t hit=0;hit<nhits;hit++){
1198         tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
1199         if (tpcHit->fQ == 0.) continue; //information about track (I.Belikov)
1200         sector=tpcHit->fSector; // sector number
1201         if(sector != isec) continue; //terminate iteration
1202         ipart=tpcHit->fTrack;
1203         particle=(TParticle*)Particles->UncheckedAt(ipart);
1204         pl=particle->Pz();
1205         pt=particle->Pt();
1206         if(pt < 1.e-9) pt=1.e-9;
1207         tanth=pl/pt;
1208         tanth = TMath::Abs(tanth);
1209         rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY);
1210         ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
1211         
1212         //   space-point resolutions
1213         
1214         sigma_rphi=SigmaY2(rpad,tanth,pt);
1215         sigma_z   =SigmaZ2(rpad,tanth   );
1216         
1217         //   cluster widths
1218         
1219         cl_rphi=ac_rphi-bc_rphi*rpad*tanth+cc_rphi*ratio*ratio;
1220         cl_z=ac_z-bc_z*rpad*tanth+cc_z*tanth*tanth;
1221         
1222         // temporary protection
1223         
1224         if(sigma_rphi < 0.) sigma_rphi=0.4e-3;
1225         if(sigma_z < 0.) sigma_z=0.4e-3;
1226         if(cl_rphi < 0.) cl_rphi=2.5e-3;
1227         if(cl_z < 0.) cl_z=2.5e-5;
1228         
1229         //
1230         
1231         //
1232         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
1233         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
1234         //
1235         //Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph;
1236         Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph;
1237         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigma_rphi));   // y
1238         xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigma_z)); // z 
1239         xyz[2]=tpcHit->fQ;                                     // q
1240         xyz[3]=sigma_rphi;                                     // fSigmaY2
1241         xyz[4]=sigma_z;                                        // fSigmaZ2
1242         
1243         Int_t tracks[5]={tpcHit->fTrack, -1, -1, sector, tpcHit->fPadRow};
1244         AddCluster(xyz,tracks);
1245         
1246       } // end of loop over hits
1247     }   // end of loop over tracks     
1248     
1249   } // end of loop over sectors  
1250   
1251 } // end of function
1252
1253 //_________________________________________________________________
1254 void AliTPC::Hits2ExactClustersSector(Int_t isec)
1255 {
1256   //--------------------------------------------------------
1257   //calculate exact cross point of track and given pad row
1258   //resulting values are expressed in "digit" coordinata
1259   //--------------------------------------------------------
1260
1261   //-----------------------------------------------------------------
1262   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
1263   //-----------------------------------------------------------------
1264   //
1265   if (fClustersArray==0){    
1266     return;
1267   }
1268   //
1269   TParticle *particle; // pointer to a given particle
1270   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1271   TClonesArray *Particles; //pointer to the particle list
1272   Int_t sector,nhits;
1273   Int_t ipart;
1274   const Int_t cmaxhits=30000;
1275   TVector * xxxx = new TVector(cmaxhits*4);
1276   TVector & xxx = *xxxx;
1277   Int_t maxhits = cmaxhits;
1278   //construct array for each padrow
1279   for (Int_t i=0; i<fTPCParam->GetNRow(isec);i++) 
1280     fClustersArray->CreateRow(isec,i);
1281   
1282   //---------------------------------------------------------------
1283   //  Get the access to the tracks 
1284   //---------------------------------------------------------------
1285   
1286   TTree *TH = gAlice->TreeH();
1287   Stat_t ntracks = TH->GetEntries();
1288   Particles=gAlice->Particles();
1289   Int_t npart = Particles->GetEntriesFast();
1290     
1291   //------------------------------------------------------------
1292   // Loop over tracks
1293   //------------------------------------------------------------
1294   
1295   for(Int_t track=0;track<ntracks;track++){
1296     ResetHits();
1297     TH->GetEvent(track);
1298     //
1299     //  Get number of the TPC hits and a pointer
1300     //  to the particles
1301     //
1302     nhits=fHits->GetEntriesFast();
1303     //
1304     // Loop over hits
1305     //
1306     Int_t currentIndex=0;
1307     Int_t lastrow=-1;  //last writen row
1308     for(Int_t hit=0;hit<nhits;hit++){
1309       tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
1310       if (tpcHit==0) continue;
1311       sector=tpcHit->fSector; // sector number
1312       if(sector != isec) continue; 
1313       ipart=tpcHit->fTrack;
1314       if (ipart<npart) particle=(TParticle*)Particles->UncheckedAt(ipart);
1315       
1316       //find row number
1317
1318       Float_t  x[3]={tpcHit->fX,tpcHit->fY,tpcHit->fZ};
1319       Int_t    index[3]={1,isec,0};
1320       Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;     
1321       if (currentrow<0) continue;
1322       if (lastrow<0) lastrow=currentrow;
1323       if (currentrow==lastrow){
1324         if ( currentIndex>=maxhits){
1325           maxhits+=cmaxhits;
1326           xxx.ResizeTo(4*maxhits);
1327         }     
1328         xxx(currentIndex*4)=x[0];
1329         xxx(currentIndex*4+1)=x[1];
1330         xxx(currentIndex*4+2)=x[2];     
1331         xxx(currentIndex*4+3)=tpcHit->fQ;
1332         currentIndex++; 
1333       }
1334       else 
1335         if (currentIndex>2){
1336           Float_t sumx=0;
1337           Float_t sumx2=0;
1338           Float_t sumx3=0;
1339           Float_t sumx4=0;
1340           Float_t sumy=0;
1341           Float_t sumxy=0;
1342           Float_t sumx2y=0;
1343           Float_t sumz=0;
1344           Float_t sumxz=0;
1345           Float_t sumx2z=0;
1346           Float_t sumq=0;
1347           for (Int_t index=0;index<currentIndex;index++){
1348             Float_t x,x2,x3,x4;
1349             x=x2=x3=x4=xxx(index*4);
1350             x2*=x;
1351             x3*=x2;
1352             x4*=x3;
1353             sumx+=x;
1354             sumx2+=x2;
1355             sumx3+=x3;
1356             sumx4+=x4;
1357             sumy+=xxx(index*4+1);
1358             sumxy+=xxx(index*4+1)*x;
1359             sumx2y+=xxx(index*4+1)*x2;
1360             sumz+=xxx(index*4+2);
1361             sumxz+=xxx(index*4+2)*x;
1362             sumx2z+=xxx(index*4+2)*x2;   
1363             sumq+=xxx(index*4+3);
1364           }
1365           Float_t CentralPad = (fTPCParam->GetNPads(isec,lastrow)-1)/2;
1366           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
1367             sumx2*(sumx*sumx3-sumx2*sumx2);
1368           
1369           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
1370             sumx2*(sumxy*sumx3-sumx2y*sumx2);
1371           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
1372             sumx2*(sumxz*sumx3-sumx2z*sumx2);
1373           
1374           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
1375             sumx2*(sumx*sumx2y-sumx2*sumxy);
1376           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
1377             sumx2*(sumx*sumx2z-sumx2*sumxz);
1378           
1379           Float_t y=detay/det+CentralPad;
1380           Float_t z=detaz/det;  
1381           Float_t by=detby/det; //y angle
1382           Float_t bz=detbz/det; //z angle
1383           sumy/=Float_t(currentIndex);
1384           sumz/=Float_t(currentIndex);
1385           AliCluster cl;
1386           cl.fX=z;
1387           cl.fY=y;
1388           cl.fQ=sumq;
1389           cl.fSigmaX2=bz;
1390           cl.fSigmaY2=by;
1391           cl.fTracks[0]=ipart;
1392           
1393           AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow));
1394           if (row!=0) row->InsertCluster(&cl);
1395           currentIndex=0;
1396           lastrow=currentrow;
1397         } //end of calculating cluster for given row
1398         
1399         
1400         
1401     } // end of loop over hits
1402   }   // end of loop over tracks 
1403   //write padrows to tree 
1404   for (Int_t ii=0; ii<fTPCParam->GetNRow(isec);ii++) {
1405     fClustersArray->StoreRow(isec,ii);    
1406     fClustersArray->ClearRow(isec,ii);        
1407   }
1408   xxxx->Delete();
1409  
1410 }
1411
1412 //__________________________________________________________________  
1413 void AliTPC::Hits2Digits()  
1414
1415  //----------------------------------------------------
1416  // Loop over all sectors
1417  //----------------------------------------------------
1418
1419   if(fTPCParam == 0){
1420     printf("AliTPCParam MUST be created firstly\n");
1421     return;
1422   } 
1423
1424  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
1425
1426 }
1427
1428
1429 //_____________________________________________________________________________
1430 void AliTPC::Hits2DigitsSector(Int_t isec)
1431 {
1432   //-------------------------------------------------------------------
1433   // TPC conversion from hits to digits.
1434   //------------------------------------------------------------------- 
1435
1436   //-----------------------------------------------------------------
1437   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1438   //-----------------------------------------------------------------
1439
1440   //-------------------------------------------------------
1441   //  Get the access to the track hits
1442   //-------------------------------------------------------
1443
1444
1445   TTree *TH = gAlice->TreeH(); // pointer to the hits tree
1446   Stat_t ntracks = TH->GetEntries();
1447
1448   if( ntracks > 0){
1449
1450   //------------------------------------------- 
1451   //  Only if there are any tracks...
1452   //-------------------------------------------
1453
1454     TObjArray **row;
1455     
1456       printf("*** Processing sector number %d ***\n",isec);
1457
1458       Int_t nrows =fTPCParam->GetNRow(isec);
1459
1460       row= new TObjArray* [nrows];
1461     
1462       MakeSector(isec,nrows,TH,ntracks,row);
1463
1464       //--------------------------------------------------------
1465       //   Digitize this sector, row by row
1466       //   row[i] is the pointer to the TObjArray of TVectors,
1467       //   each one containing electrons accepted on this
1468       //   row, assigned into tracks
1469       //--------------------------------------------------------
1470
1471       Int_t i;
1472
1473       if (fDigitsArray->GetTree()==0) fDigitsArray->MakeTree();
1474
1475       for (i=0;i<nrows;i++){
1476
1477         AliDigits * dig = fDigitsArray->CreateRow(isec,i); 
1478
1479         DigitizeRow(i,isec,row);
1480
1481         fDigitsArray->StoreRow(isec,i);
1482
1483         Int_t ndig = dig->GetSize(); 
1484  
1485         printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);
1486         
1487         fDigitsArray->ClearRow(isec,i);  
1488
1489    
1490        } // end of the sector digitization
1491
1492       for(i=0;i<nrows;i++){
1493         row[i]->Delete();     
1494       }
1495       
1496        delete [] row; // delete the array of pointers to TObjArray-s
1497         
1498   } // ntracks >0
1499
1500 } // end of Hits2DigitsSector
1501
1502
1503 //_____________________________________________________________________________
1504 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1505 {
1506   //-----------------------------------------------------------
1507   // Single row digitization, coupling from the neighbouring
1508   // rows taken into account
1509   //-----------------------------------------------------------
1510
1511   //-----------------------------------------------------------------
1512   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1513   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1514   //-----------------------------------------------------------------
1515  
1516
1517   Float_t zerosup = fTPCParam->GetZeroSup();
1518   Int_t nrows =fTPCParam->GetNRow(isec);
1519   fCurrentIndex[1]= isec;
1520   
1521
1522   Int_t n_of_pads = fTPCParam->GetNPads(isec,irow);
1523   Int_t n_of_tbins = fTPCParam->GetMaxTBin();
1524   Int_t IndexRange[4];
1525   //
1526   //  Integrated signal for this row
1527   //  and a single track signal
1528   //    
1529   TMatrix *m1   = new TMatrix(0,n_of_pads,0,n_of_tbins); // integrated
1530   TMatrix *m2   = new TMatrix(0,n_of_pads,0,n_of_tbins); // single
1531   //
1532   TMatrix &Total  = *m1;
1533
1534   //  Array of pointers to the label-signal list
1535
1536   Int_t NofDigits = n_of_pads*n_of_tbins; // number of digits for this row
1537   Float_t  **pList = new Float_t* [NofDigits]; 
1538
1539   Int_t lp;
1540   Int_t i1;   
1541   for(lp=0;lp<NofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1542   //
1543   //calculate signal 
1544   //
1545   Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
1546   Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
1547   for (Int_t row= row1;row<=row2;row++){
1548     Int_t nTracks= rows[row]->GetEntries();
1549     for (i1=0;i1<nTracks;i1++){
1550       fCurrentIndex[2]= row;
1551       fCurrentIndex[3]=irow;
1552       if (row==irow){
1553         m2->Zero();  // clear single track signal matrix
1554         Float_t TrackLabel = GetSignal(rows[row],i1,m2,m1,IndexRange); 
1555         GetList(TrackLabel,n_of_pads,m2,IndexRange,pList);
1556       }
1557       else   GetSignal(rows[row],i1,0,m1,IndexRange);
1558     }
1559   }
1560          
1561   Int_t tracks[3];
1562
1563   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1564   for(Int_t ip=0;ip<n_of_pads;ip++){
1565     for(Int_t it=0;it<n_of_tbins;it++){
1566
1567       Float_t q = Total(ip,it);
1568
1569       Int_t gi =it*n_of_pads+ip; // global index
1570
1571       q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()); 
1572
1573       q = (Int_t)q;
1574
1575       if(q <=zerosup) continue; // do not fill zeros
1576       if(q > adc_sat) q = adc_sat;  // saturation
1577
1578       //
1579       //  "real" signal or electronic noise (list = -1)?
1580       //    
1581
1582       for(Int_t j1=0;j1<3;j1++){
1583         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1;
1584       }
1585
1586 //Begin_Html
1587 /*
1588   <A NAME="AliDigits"></A>
1589   using of AliDigits object
1590 */
1591 //End_Html
1592       dig->SetDigitFast((Short_t)q,it,ip);
1593       if (fDigitsArray->IsSimulated())
1594         {
1595          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1596          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1597          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1598         }
1599      
1600     
1601     } // end of loop over time buckets
1602   }  // end of lop over pads 
1603
1604   //
1605   //  This row has been digitized, delete nonused stuff
1606   //
1607
1608   for(lp=0;lp<NofDigits;lp++){
1609     if(pList[lp]) delete [] pList[lp];
1610   }
1611   
1612   delete [] pList;
1613
1614   delete m1;
1615   delete m2;
1616   //  delete m3;
1617
1618 } // end of DigitizeRow
1619
1620 //_____________________________________________________________________________
1621
1622 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
1623                           Int_t *IndexRange)
1624 {
1625
1626   //---------------------------------------------------------------
1627   //  Calculates 2-D signal (pad,time) for a single track,
1628   //  returns a pointer to the signal matrix and the track label 
1629   //  No digitization is performed at this level!!!
1630   //---------------------------------------------------------------
1631
1632   //-----------------------------------------------------------------
1633   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1634   // Modified: Marian Ivanov 
1635   //-----------------------------------------------------------------
1636
1637   TVector *tv;
1638  
1639   tv = (TVector*)p1->At(ntr); // pointer to a track
1640   TVector &v = *tv;
1641   
1642   Float_t label = v(0);
1643   Int_t CentralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3])-1)/2;
1644
1645   Int_t nElectrons = (tv->GetNrows()-1)/4;
1646   IndexRange[0]=9999; // min pad
1647   IndexRange[1]=-1; // max pad
1648   IndexRange[2]=9999; //min time
1649   IndexRange[3]=-1; // max time
1650
1651   //  Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
1652
1653   TMatrix &signal = *m1;
1654   TMatrix &total = *m2;
1655   //
1656   //  Loop over all electrons
1657   //
1658   for(Int_t nel=0; nel<nElectrons; nel++){
1659     Int_t idx=nel*4;
1660     Float_t aval =  v(idx+4);
1661     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac(); 
1662     Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
1663     Int_t n = fTPCParam->CalcResponse(xyz,fCurrentIndex,fCurrentIndex[3]);
1664     
1665     if (n>0) for (Int_t i =0; i<n; i++){
1666        Int_t *index = fTPCParam->GetResBin(i);        
1667        Int_t pad=index[1]+CentralPad;  //in digit coordinates central pad has coordinate 0
1668        if ( ( pad<(fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]))) && (pad>0)) {
1669          Int_t time=index[2];    
1670          Float_t weight = fTPCParam->GetResWeight(i); //we normalise response to ADC channel
1671          weight *= eltoadcfac;
1672          
1673          if (m1!=0) signal(pad,time)+=weight; 
1674          total(pad,time)+=weight;
1675          IndexRange[0]=TMath::Min(IndexRange[0],pad);
1676          IndexRange[1]=TMath::Max(IndexRange[1],pad);
1677          IndexRange[2]=TMath::Min(IndexRange[2],time);
1678          IndexRange[3]=TMath::Max(IndexRange[3],time); 
1679        }         
1680     }
1681   } // end of loop over electrons
1682   
1683   return label; // returns track label when finished
1684 }
1685
1686 //_____________________________________________________________________________
1687 void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *IndexRange,
1688                      Float_t **pList)
1689 {
1690   //----------------------------------------------------------------------
1691   //  Updates the list of tracks contributing to digits for a given row
1692   //----------------------------------------------------------------------
1693
1694   //-----------------------------------------------------------------
1695   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1696   //-----------------------------------------------------------------
1697
1698   TMatrix &signal = *m;
1699
1700   // lop over nonzero digits
1701
1702   for(Int_t it=IndexRange[2];it<IndexRange[3]+1;it++){
1703     for(Int_t ip=IndexRange[0];ip<IndexRange[1]+1;ip++){
1704
1705
1706         // accept only the contribution larger than 500 electrons (1/2 s_noise)
1707
1708         if(signal(ip,it)<0.5) continue; 
1709
1710
1711         Int_t GlobalIndex = it*np+ip; // GlobalIndex starts from 0!
1712         
1713         if(!pList[GlobalIndex]){
1714         
1715           // 
1716           // Create new list (6 elements - 3 signals and 3 labels),
1717           //
1718
1719           pList[GlobalIndex] = new Float_t [6];
1720
1721           // set list to -1 
1722
1723           *pList[GlobalIndex] = -1.;
1724           *(pList[GlobalIndex]+1) = -1.;
1725           *(pList[GlobalIndex]+2) = -1.;
1726           *(pList[GlobalIndex]+3) = -1.;
1727           *(pList[GlobalIndex]+4) = -1.;
1728           *(pList[GlobalIndex]+5) = -1.;
1729
1730
1731           *pList[GlobalIndex] = label;
1732           *(pList[GlobalIndex]+3) = signal(ip,it);
1733         }
1734         else{
1735
1736           // check the signal magnitude
1737
1738           Float_t highest = *(pList[GlobalIndex]+3);
1739           Float_t middle = *(pList[GlobalIndex]+4);
1740           Float_t lowest = *(pList[GlobalIndex]+5);
1741
1742           //
1743           //  compare the new signal with already existing list
1744           //
1745
1746           if(signal(ip,it)<lowest) continue; // neglect this track
1747
1748           //
1749
1750           if (signal(ip,it)>highest){
1751             *(pList[GlobalIndex]+5) = middle;
1752             *(pList[GlobalIndex]+4) = highest;
1753             *(pList[GlobalIndex]+3) = signal(ip,it);
1754
1755             *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1);
1756             *(pList[GlobalIndex]+1) = *pList[GlobalIndex];
1757             *pList[GlobalIndex] = label;
1758           }
1759           else if (signal(ip,it)>middle){
1760             *(pList[GlobalIndex]+5) = middle;
1761             *(pList[GlobalIndex]+4) = signal(ip,it);
1762
1763             *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1);
1764             *(pList[GlobalIndex]+1) = label;
1765           }
1766           else{
1767             *(pList[GlobalIndex]+5) = signal(ip,it);
1768             *(pList[GlobalIndex]+2) = label;
1769           }
1770         }
1771
1772     } // end of loop over pads
1773   } // end of loop over time bins
1774
1775
1776
1777 }//end of GetList
1778 //___________________________________________________________________
1779 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1780                         Stat_t ntracks,TObjArray **row)
1781 {
1782
1783   //-----------------------------------------------------------------
1784   // Prepares the sector digitization, creates the vectors of
1785   // tracks for each row of this sector. The track vector
1786   // contains the track label and the position of electrons.
1787   //-----------------------------------------------------------------
1788
1789   //-----------------------------------------------------------------
1790   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1791   //-----------------------------------------------------------------
1792
1793   Float_t gasgain = fTPCParam->GetGasGain();
1794   Int_t i;
1795   Float_t xyz[4]; 
1796
1797   AliTPChit *tpcHit; // pointer to a sigle TPC hit    
1798  
1799   //----------------------------------------------
1800   // Create TObjArray-s, one for each row,
1801   // each TObjArray will store the TVectors
1802   // of electrons, one TVector per each track.
1803   //---------------------------------------------- 
1804     
1805   for(i=0; i<nrows; i++){
1806     row[i] = new TObjArray;
1807   }
1808   Int_t *n_of_electrons = new Int_t [nrows]; // electron counter for each row
1809   TVector **tracks = new TVector* [nrows]; //pointers to the track vectors
1810
1811   //--------------------------------------------------------------------
1812   //  Loop over tracks, the "track" contains the full history
1813   //--------------------------------------------------------------------
1814
1815   Int_t previousTrack,currentTrack;
1816   previousTrack = -1; // nothing to store so far!
1817
1818   for(Int_t track=0;track<ntracks;track++){
1819
1820     ResetHits();
1821
1822     TH->GetEvent(track); // get next track
1823     Int_t nhits = fHits->GetEntriesFast(); // get number of hits for this track
1824
1825     if(nhits == 0) continue; // no hits in the TPC for this track
1826
1827     //--------------------------------------------------------------
1828     //  Loop over hits
1829     //--------------------------------------------------------------
1830
1831     for(Int_t hit=0;hit<nhits;hit++){
1832
1833       tpcHit = (AliTPChit*)fHits->UncheckedAt(hit); // get a pointer to a hit
1834       
1835       Int_t sector=tpcHit->fSector; // sector number
1836       if(sector != isec) continue; 
1837
1838         currentTrack = tpcHit->fTrack; // track number
1839         if(currentTrack != previousTrack){
1840                           
1841            // store already filled fTrack
1842               
1843            for(i=0;i<nrows;i++){
1844              if(previousTrack != -1){
1845                if(n_of_electrons[i]>0){
1846                  TVector &v = *tracks[i];
1847                  v(0) = previousTrack;
1848                  tracks[i]->ResizeTo(4*n_of_electrons[i]+1); // shrink if necessary
1849                  row[i]->Add(tracks[i]);                     
1850                }
1851                else{
1852                  delete tracks[i]; // delete empty TVector
1853                  tracks[i]=0;
1854                }
1855              }
1856
1857              n_of_electrons[i]=0;
1858              tracks[i] = new TVector(481); // TVectors for the next fTrack
1859
1860            } // end of loop over rows
1861                
1862            previousTrack=currentTrack; // update track label 
1863         }
1864            
1865         Int_t QI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
1866
1867        //---------------------------------------------------
1868        //  Calculate the electron attachment probability
1869        //---------------------------------------------------
1870
1871
1872         Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->fZ))
1873                                                         /fTPCParam->GetDriftV(); 
1874         // in microseconds!     
1875         Float_t AttProb = fTPCParam->GetAttCoef()*
1876           fTPCParam->GetOxyCont()*time; //  fraction! 
1877    
1878         //-----------------------------------------------
1879         //  Loop over electrons
1880         //-----------------------------------------------
1881         Int_t index[3];
1882         index[1]=isec;
1883         for(Int_t nel=0;nel<QI;nel++){
1884           // skip if electron lost due to the attachment
1885           if((gRandom->Rndm(0)) < AttProb) continue; // electron lost!
1886           xyz[0]=tpcHit->fX;
1887           xyz[1]=tpcHit->fY;
1888           xyz[2]=tpcHit->fZ;      
1889           xyz[3]= (Float_t) (-gasgain*TMath::Log(gRandom->Rndm()));
1890           index[0]=1;
1891           
1892           TransportElectron(xyz,index); //MI change -august       
1893           Int_t row_number;
1894           fTPCParam->GetPadRow(xyz,index); //MI change august
1895           row_number = index[2];
1896           //transform position to local digit coordinates
1897           //relative to nearest pad row 
1898           if ((row_number<0)||row_number>=fTPCParam->GetNRow(isec)) continue;     
1899           n_of_electrons[row_number]++;   
1900           //----------------------------------
1901           // Expand vector if necessary
1902           //----------------------------------
1903           if(n_of_electrons[row_number]>120){
1904             Int_t range = tracks[row_number]->GetNrows();
1905             if((n_of_electrons[row_number])>(range-1)/4){
1906         
1907               tracks[row_number]->ResizeTo(range+400); // Add 100 electrons
1908             }
1909           }
1910           
1911           TVector &v = *tracks[row_number];
1912           Int_t idx = 4*n_of_electrons[row_number]-3;
1913
1914           v(idx)=  xyz[0];   // X - pad row coordinate
1915           v(idx+1)=xyz[1];   // Y - pad coordinate (along the pad-row)
1916           v(idx+2)=xyz[2];   // Z - time bin coordinate
1917           v(idx+3)=xyz[3];   // avalanche size  
1918         } // end of loop over electrons
1919         
1920       } // end of loop over hits
1921     } // end of loop over tracks
1922
1923     //
1924     //   store remaining track (the last one) if not empty
1925     //
1926
1927      for(i=0;i<nrows;i++){
1928        if(n_of_electrons[i]>0){
1929           TVector &v = *tracks[i];
1930           v(0) = previousTrack;
1931           tracks[i]->ResizeTo(4*n_of_electrons[i]+1); // shrink if necessary
1932           row[i]->Add(tracks[i]);  
1933         }
1934         else{
1935           delete tracks[i];
1936           tracks[i]=0;
1937         }  
1938       }  
1939
1940           delete [] tracks;
1941           delete [] n_of_electrons;
1942  
1943
1944 } // end of MakeSector
1945
1946
1947 //_____________________________________________________________________________
1948 void AliTPC::Init()
1949 {
1950   //
1951   // Initialise TPC detector after definition of geometry
1952   //
1953   Int_t i;
1954   //
1955   printf("\n");
1956   for(i=0;i<35;i++) printf("*");
1957   printf(" TPC_INIT ");
1958   for(i=0;i<35;i++) printf("*");
1959   printf("\n");
1960   //
1961   for(i=0;i<80;i++) printf("*");
1962   printf("\n");
1963 }
1964
1965 //_____________________________________________________________________________
1966 void AliTPC::MakeBranch(Option_t* option)
1967 {
1968   //
1969   // Create Tree branches for the TPC.
1970   //
1971   Int_t buffersize = 4000;
1972   char branchname[10];
1973   sprintf(branchname,"%s",GetName());
1974
1975   AliDetector::MakeBranch(option);
1976
1977   char *D = strstr(option,"D");
1978
1979   if (fDigits   && gAlice->TreeD() && D) {
1980     gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
1981     printf("Making Branch %s for digits\n",branchname);
1982   }     
1983
1984   char *R = strstr(option,"R");
1985
1986   if (fClusters && gAlice->TreeR() && R) {
1987     gAlice->TreeR()->Branch(branchname,&fClusters, buffersize);
1988     printf("Making Branch %s for Clusters\n",branchname);
1989   }     
1990 }
1991  
1992 //_____________________________________________________________________________
1993 void AliTPC::ResetDigits()
1994 {
1995   //
1996   // Reset number of digits and the digits array for this detector
1997   // reset clusters
1998   //
1999   fNdigits   = 0;
2000   if (fDigits)   fDigits->Clear();
2001   fNclusters = 0;
2002   if (fClusters) fClusters->Clear();
2003 }
2004
2005 //_____________________________________________________________________________
2006 void AliTPC::SetSecAL(Int_t sec)
2007 {
2008   //---------------------------------------------------
2009   // Activate/deactivate selection for lower sectors
2010   //---------------------------------------------------
2011
2012   //-----------------------------------------------------------------
2013   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2014   //-----------------------------------------------------------------
2015
2016   fSecAL = sec;
2017 }
2018
2019 //_____________________________________________________________________________
2020 void AliTPC::SetSecAU(Int_t sec)
2021 {
2022   //----------------------------------------------------
2023   // Activate/deactivate selection for upper sectors
2024   //---------------------------------------------------
2025
2026   //-----------------------------------------------------------------
2027   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2028   //-----------------------------------------------------------------
2029
2030   fSecAU = sec;
2031 }
2032
2033 //_____________________________________________________________________________
2034 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
2035 {
2036   //----------------------------------------
2037   // Select active lower sectors
2038   //----------------------------------------
2039
2040   //-----------------------------------------------------------------
2041   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2042   //-----------------------------------------------------------------
2043
2044   fSecLows[0] = s1;
2045   fSecLows[1] = s2;
2046   fSecLows[2] = s3;
2047   fSecLows[3] = s4;
2048   fSecLows[4] = s5;
2049   fSecLows[5] = s6;
2050 }
2051
2052 //_____________________________________________________________________________
2053 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
2054                        Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10, 
2055                        Int_t s11 , Int_t s12)
2056 {
2057   //--------------------------------
2058   // Select active upper sectors
2059   //--------------------------------
2060
2061   //-----------------------------------------------------------------
2062   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2063   //-----------------------------------------------------------------
2064
2065   fSecUps[0] = s1;
2066   fSecUps[1] = s2;
2067   fSecUps[2] = s3;
2068   fSecUps[3] = s4;
2069   fSecUps[4] = s5;
2070   fSecUps[5] = s6;
2071   fSecUps[6] = s7;
2072   fSecUps[7] = s8;
2073   fSecUps[8] = s9;
2074   fSecUps[9] = s10;
2075   fSecUps[10] = s11;
2076   fSecUps[11] = s12;
2077 }
2078
2079 //_____________________________________________________________________________
2080 void AliTPC::SetSens(Int_t sens)
2081 {
2082
2083   //-------------------------------------------------------------
2084   // Activates/deactivates the sensitive strips at the center of
2085   // the pad row -- this is for the space-point resolution calculations
2086   //-------------------------------------------------------------
2087
2088   //-----------------------------------------------------------------
2089   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2090   //-----------------------------------------------------------------
2091
2092   fSens = sens;
2093 }
2094  
2095 void AliTPC::SetSide(Float_t side)
2096 {
2097   fSide = side;
2098  
2099 }
2100 //____________________________________________________________________________
2101 void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
2102                            Float_t p2,Float_t p3)
2103 {
2104
2105  fNoComp = nc;
2106  
2107  fMixtComp[0]=c1;
2108  fMixtComp[1]=c2;
2109  fMixtComp[2]=c3;
2110
2111  fMixtProp[0]=p1;
2112  fMixtProp[1]=p2;
2113  fMixtProp[2]=p3; 
2114  
2115  
2116 }
2117 //_____________________________________________________________________________
2118
2119 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2120 {
2121   //
2122   // electron transport taking into account:
2123   // 1. diffusion, 
2124   // 2.ExB at the wires
2125   // 3. nonisochronity
2126   //
2127   // xyz and index must be already transformed to system 1
2128   //
2129
2130   fTPCParam->Transform1to2(xyz,index);
2131   
2132   //add diffusion
2133   Float_t driftl=xyz[2];
2134   if(driftl<0.01) driftl=0.01;
2135   driftl=TMath::Sqrt(driftl);
2136   Float_t sig_t = driftl*(fTPCParam->GetDiffT());
2137   Float_t sig_l = driftl*(fTPCParam->GetDiffL());
2138   xyz[0]=gRandom->Gaus(xyz[0],sig_t);
2139   xyz[1]=gRandom->Gaus(xyz[1],sig_t);
2140   xyz[2]=gRandom->Gaus(xyz[2],sig_l);
2141
2142   // ExB
2143   
2144   if (fTPCParam->GetMWPCReadout()==kTRUE){
2145     Float_t x1=xyz[0];
2146     fTPCParam->Transform2to2NearestWire(xyz,index);
2147     Float_t dx=xyz[0]-x1;
2148     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2149   }
2150   //add nonisochronity (not implemented yet)
2151   
2152 }
2153 //_____________________________________________________________________________
2154 void AliTPC::Streamer(TBuffer &R__b)
2155 {
2156   //
2157   // Stream an object of class AliTPC.
2158   //
2159    if (R__b.IsReading()) {
2160       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
2161       AliDetector::Streamer(R__b);
2162       if (R__v < 2) return;
2163       R__b >> fNsectors;
2164       R__b >> fNclusters;
2165       R__b >> fNtracks;
2166
2167    } else {
2168       R__b.WriteVersion(AliTPC::IsA());
2169       AliDetector::Streamer(R__b);
2170       R__b << fNsectors;
2171       R__b << fNclusters;
2172       R__b << fNtracks;
2173    }
2174 }
2175   
2176 ClassImp(AliTPCcluster)
2177  
2178 //_____________________________________________________________________________
2179 AliTPCcluster::AliTPCcluster(Float_t *hits, Int_t *lab)
2180 {
2181   //
2182   // Creates a simulated cluster for the TPC
2183   //
2184   fTracks[0]  = lab[0];
2185   fTracks[1]  = lab[1];
2186   fTracks[2]  = lab[2];
2187   fSector     = lab[3];
2188   fPadRow     = lab[4];
2189   fY          = hits[0];
2190   fZ          = hits[1];
2191   fQ          = hits[2];
2192   fSigmaY2    = hits[3];
2193   fSigmaZ2    = hits[4];
2194 }
2195  
2196 //_____________________________________________________________________________
2197 void AliTPCcluster::GetXYZ(Float_t *x, const AliTPCParam *par) const 
2198 {
2199   //
2200   // Transformation from local to global coordinate system
2201   //
2202   x[0]=par->GetPadRowRadii(fSector,fPadRow);
2203   x[1]=fY;
2204   x[2]=fZ;
2205   Float_t cs, sn, tmp;
2206   par->AdjustCosSin(fSector,cs,sn);
2207   tmp = x[0]*cs-x[1]*sn;
2208   x[1]= x[0]*sn+x[1]*cs; x[0]=tmp;
2209 }
2210  
2211 //_____________________________________________________________________________
2212 Int_t AliTPCcluster::Compare(TObject * o)
2213 {
2214   //
2215   // compare two clusters according y coordinata
2216   //
2217   AliTPCcluster *cl= (AliTPCcluster *)o;
2218   if (fY<cl->fY) return -1;
2219   if (fY==cl->fY) return 0;
2220   return 1;  
2221 }
2222
2223 Bool_t AliTPCcluster::IsSortable() const
2224 {
2225   //
2226   //make AliTPCcluster sortabale
2227   //
2228   return kTRUE; 
2229 }
2230
2231
2232
2233 ClassImp(AliTPCdigit)
2234  
2235 //_____________________________________________________________________________
2236 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
2237   AliDigit(tracks)
2238 {
2239   //
2240   // Creates a TPC digit object
2241   //
2242   fSector     = digits[0];
2243   fPadRow     = digits[1];
2244   fPad        = digits[2];
2245   fTime       = digits[3];
2246   fSignal     = digits[4];
2247 }
2248
2249  
2250 ClassImp(AliTPChit)
2251  
2252 //_____________________________________________________________________________
2253 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
2254 AliHit(shunt,track)
2255 {
2256   //
2257   // Creates a TPC hit object
2258   //
2259   fSector     = vol[0];
2260   fPadRow     = vol[1];
2261   fX          = hits[0];
2262   fY          = hits[1];
2263   fZ          = hits[2];
2264   fQ          = hits[3];
2265 }
2266  
2267  
2268 ClassImp(AliTPCtrack)
2269  
2270 //_____________________________________________________________________________
2271 AliTPCtrack::AliTPCtrack(Float_t *hits)
2272 {
2273   //
2274   // Default creator for a TPC reconstructed track object
2275   //
2276   fX=hits[0]; // This is dummy code !
2277 }
2278 //_________________________________________________________________________
2279
2280 AliTPCtrack::AliTPCtrack(const AliTPCcluster *c,const TVector& xx,
2281                          const TMatrix& CC, Double_t xref, Double_t alpha):
2282   x(xx),C(CC),fClusters(200)
2283 {
2284   //-----------------------------------------------------------------
2285   // This is the main track constructor.
2286   //
2287   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2288   //-----------------------------------------------------------------
2289   fX=xref;
2290   fAlpha=alpha;
2291   fChi2=0.;
2292   fClusters.AddLast((AliTPCcluster*)(c));
2293 }
2294
2295 //_____________________________________________________________________________
2296 AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : x(t.x), C(t.C),
2297   fClusters(t.fClusters.GetEntriesFast()) 
2298 {
2299   //-----------------------------------------------------------------
2300   // This is a track copy constructor.
2301   //
2302   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2303   //-----------------------------------------------------------------
2304   fX=t.fX;
2305   fChi2=t.fChi2;
2306   fAlpha=t.fAlpha;
2307   Int_t n=t.fClusters.GetEntriesFast();
2308   for (Int_t i=0; i<n; i++) fClusters.AddLast(t.fClusters.UncheckedAt(i));
2309 }
2310
2311 //_____________________________________________________________________________
2312 Int_t AliTPCtrack::Compare(TObject *o) {
2313   //-----------------------------------------------------------------
2314   // This function compares tracks according to the uncertainty of their
2315   // position in Y.
2316   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2317   //-----------------------------------------------------------------
2318   AliTPCtrack *t=(AliTPCtrack*)o;
2319   Double_t co=t->GetSigmaY2();
2320   Double_t c =GetSigmaY2();
2321   if (c>co) return 1;
2322   else if (c<co) return -1;
2323   return 0;
2324 }
2325
2326 //_____________________________________________________________________________
2327 Int_t AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
2328 {
2329   //-----------------------------------------------------------------
2330   // This function propagates a track to a reference plane x=xk.
2331   //
2332   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2333   //-----------------------------------------------------------------
2334   if (TMath::Abs(x(2)*xk - x(3)) >= 0.999) {
2335     if (*this>4) cerr<<*this<<" AliTPCtrack warning: Propagation failed !\n";
2336     return 0;
2337   }
2338
2339   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1, y1=x(0), z1=x(1);
2340   Double_t c1=x(2)*x1 - x(3), r1=sqrt(1.- c1*c1);
2341   Double_t c2=x(2)*x2 - x(3), r2=sqrt(1.- c2*c2);
2342   
2343   x(0) += dx*(c1+c2)/(r1+r2);
2344   x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
2345
2346   TMatrix F(5,5); F.UnitMatrix();
2347   Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2;
2348   F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
2349   F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
2350   Double_t cr=c1*r2+c2*r1;
2351   F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
2352   F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
2353   F(1,4)= dx*cc/cr; 
2354   TMatrix tmp(F,TMatrix::kMult,C);
2355   C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2356
2357   fX=x2;
2358
2359   //Multiple scattering******************
2360   Double_t ey=x(2)*fX - x(3);
2361   Double_t ex=sqrt(1-ey*ey);
2362   Double_t ez=x(4);
2363   TMatrix Q(5,5); Q=0.;
2364   Q(2,2)=ez*ez+ey*ey;   Q(2,3)=-ex*ey;       Q(2,4)=-ex*ez;
2365   Q(3,2)=Q(2,3);        Q(3,3)= ez*ez+ex*ex; Q(3,4)=-ey*ez;
2366   Q(4,2)=Q(2,4);        Q(4,3)= Q(3,4);      Q(4,4)=1.;
2367   
2368   F=0;
2369   F(2,2)=-x(2)*ex;          F(2,3)=-x(2)*ey;
2370   F(3,2)=-ex*(x(2)*fX-ey);  F(3,3)=-(1.+ x(2)*fX*ey - ey*ey);
2371   F(4,2)=-ez*ex;            F(4,3)=-ez*ey;           F(4,4)=1.;
2372   
2373   tmp.Mult(F,Q);
2374   Q.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2375   
2376   Double_t p2=GetPt()*GetPt()*(1.+x(4)*x(4));
2377   Double_t beta2=p2/(p2 + pm*pm);
2378   Double_t d=sqrt((x1-fX)*(x1-fX)+(y1-x(0))*(y1-x(0))+(z1-x(1))*(z1-x(1)));
2379   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
2380   Q*=theta2;
2381   C+=Q;
2382
2383   //Energy losses************************
2384   Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho;
2385   if (x1 < x2) dE=-dE;
2386   x(2)*=(1.- sqrt(p2+pm*pm)/p2*dE);
2387   //x(3)*=(1.- sqrt(p2+pm*pm)/p2*dE);
2388
2389   return 1;
2390 }
2391
2392 //_____________________________________________________________________________
2393 void AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm) 
2394 {
2395   //-----------------------------------------------------------------
2396   // This function propagates tracks to the "vertex".
2397   //
2398   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2399   //-----------------------------------------------------------------
2400   Double_t c=x(2)*fX - x(3);
2401   Double_t tgf=-x(3)/(x(2)*x(0) + sqrt(1-c*c));
2402   Double_t snf=tgf/sqrt(1.+ tgf*tgf);
2403   Double_t xv=(x(3)+snf)/x(2);
2404   PropagateTo(xv,x0,rho,pm);
2405 }
2406
2407 //_____________________________________________________________________________
2408 void AliTPCtrack::Update(const AliTPCcluster *c, Double_t chisq)
2409 {
2410   //-----------------------------------------------------------------
2411   // This function associates a clusters with this track.
2412   //
2413   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2414   //-----------------------------------------------------------------
2415   TMatrix H(2,5); H.UnitMatrix();
2416   TMatrix Ht(TMatrix::kTransposed,H);
2417   TVector m(2);   m(0)=c->fY; m(1)=c->fZ; 
2418   TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
2419
2420   TMatrix tmp(H,TMatrix::kMult,C);
2421   TMatrix R(tmp,TMatrix::kMult,Ht); R+=V;
2422   
2423   Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
2424   R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1); 
2425   R(1,0)*=-1; R(0,1)=R(1,0);
2426   R*=1./det;
2427   
2428   //R.Invert();
2429   
2430   TMatrix K(C,TMatrix::kMult,Ht); K*=R;
2431   
2432   TVector savex=x;
2433   x*=H; x-=m;
2434
2435   x*=-1; x*=K; x+=savex;
2436   if (TMath::Abs(x(2)*fX-x(3)) >= 0.999) {
2437     if (*this>4) cerr<<*this<<" AliTPCtrack warning: Filtering failed !\n";
2438     x=savex;
2439     return;
2440   }
2441   
2442   TMatrix saveC=C;
2443   C.Mult(K,tmp); C-=saveC; C*=-1;
2444
2445   fClusters.AddLast((AliTPCcluster*)c);
2446   fChi2 += chisq;
2447 }
2448
2449 //_____________________________________________________________________________
2450 Int_t AliTPCtrack::Rotate(Double_t alpha)
2451 {
2452   //-----------------------------------------------------------------
2453   // This function rotates this track.
2454   //
2455   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2456   //-----------------------------------------------------------------
2457   fAlpha += alpha;
2458   
2459   Double_t x1=fX, y1=x(0);
2460   Double_t ca=cos(alpha), sa=sin(alpha);
2461   Double_t r1=x(2)*fX - x(3);
2462   
2463   fX = x1*ca + y1*sa;
2464   x(0)=-x1*sa + y1*ca;
2465   x(3)=x(3)*ca + (x(2)*y1 + sqrt(1.- r1*r1))*sa;
2466   
2467   Double_t r2=x(2)*fX - x(3);
2468   if (TMath::Abs(r2) >= 0.999) {
2469     if (*this>4) cerr<<*this<<" AliTPCtrack warning: Rotation failed !\n";
2470     return 0;
2471   }
2472   
2473   Double_t y0=x(0) + sqrt(1.- r2*r2)/x(2);
2474   if ((x(0)-y0)*x(2) >= 0.) {
2475     if (*this>4) cerr<<*this<<" AliTPCtrack warning: Rotation failed !!!\n";
2476     return 0;
2477   }
2478   
2479   TMatrix F(5,5); F.UnitMatrix();
2480   F(0,0)=ca;
2481   F(3,0)=x(2)*sa; 
2482   F(3,2)=(y1 - r1*x1/sqrt(1.- r1*r1))*sa; 
2483   F(3,3)= ca + sa*r1/sqrt(1.- r1*r1);
2484   TMatrix tmp(F,TMatrix::kMult,C); 
2485   // Double_t dy2=C(0,0);
2486   C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2487   // C(0,0)+=dy2*sa*sa*r1*r1/(1.- r1*r1);
2488   // C(1,1)+=dy2*sa*sa*x(4)*x(4)/(1.- r1*r1);
2489   
2490   return 1;
2491 }
2492
2493 //_____________________________________________________________________________
2494 void AliTPCtrack::UseClusters() const 
2495 {
2496   //-----------------------------------------------------------------
2497   // This function marks clusters associated with this track.
2498   //
2499   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2500   //-----------------------------------------------------------------
2501   Int_t num_of_clusters=fClusters.GetEntriesFast();
2502   for (Int_t i=0; i<num_of_clusters; i++) {
2503     //if (i<=14) continue;
2504     AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(i);
2505     c->Use();   
2506   }
2507 }
2508
2509 //_____________________________________________________________________________
2510 Double_t AliTPCtrack::GetPredictedChi2(const AliTPCcluster *c) const 
2511 {
2512   //-----------------------------------------------------------------
2513   // This function calculates a predicted chi2 increment.
2514   //
2515   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2516   //-----------------------------------------------------------------
2517   TMatrix H(2,5); H.UnitMatrix();
2518   TVector m(2);   m(0)=c->fY; m(1)=c->fZ; 
2519   TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
2520   TVector res=x;  res*=H; res-=m; //res*=-1; 
2521   TMatrix tmp(H,TMatrix::kMult,C);
2522   TMatrix R(tmp,TMatrix::kMult,TMatrix(TMatrix::kTransposed,H)); R+=V;
2523   
2524   Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
2525   if (TMath::Abs(det) < 1.e-10) {
2526     if (*this>4) cerr<<*this<<" AliTPCtrack warning: Singular matrix !\n";
2527     return 1e10;
2528   }
2529   R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1); 
2530   R(1,0)*=-1; R(0,1)=R(1,0);
2531   R*=1./det;
2532   
2533   //R.Invert();
2534   
2535   TVector r=res;
2536   res*=R;
2537   return r*res;
2538 }
2539
2540 //_____________________________________________________________________________
2541 struct S { Int_t lab; Int_t max; };
2542 Int_t AliTPCtrack::GetLabel(Int_t nrows) const 
2543 {
2544   //-----------------------------------------------------------------
2545   // This function returns the track label. If label<0, this track is fake.
2546   //
2547   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2548   //-----------------------------------------------------------------
2549   Int_t num_of_clusters=fClusters.GetEntriesFast();
2550   S *s=new S[num_of_clusters];
2551   Int_t i;
2552   for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
2553   
2554   Int_t lab=123456789;
2555   for (i=0; i<num_of_clusters; i++) {
2556     AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(i);
2557     lab=TMath::Abs(c->fTracks[0]);
2558     Int_t j;
2559     for (j=0; j<num_of_clusters; j++)
2560       if (s[j].lab==lab || s[j].max==0) break;
2561     s[j].lab=lab;
2562     s[j].max++;
2563   }
2564   
2565   Int_t max=0;
2566   for (i=0; i<num_of_clusters; i++) 
2567     if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
2568     
2569   delete[] s;
2570   
2571   for (i=0; i<num_of_clusters; i++) {
2572     AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(i);
2573     if (TMath::Abs(c->fTracks[1]) == lab ||
2574         TMath::Abs(c->fTracks[2]) == lab ) max++;
2575   }
2576   
2577   if (1.-Float_t(max)/num_of_clusters > 0.10) return -lab;
2578   
2579   Int_t tail=Int_t(0.08*nrows);
2580   if (num_of_clusters < tail) return lab;
2581   
2582   max=0;
2583   for (i=1; i<=tail; i++) {
2584     AliTPCcluster *c=(AliTPCcluster*)fClusters.UncheckedAt(num_of_clusters-i);
2585     if (lab == TMath::Abs(c->fTracks[0]) ||
2586         lab == TMath::Abs(c->fTracks[1]) ||
2587         lab == TMath::Abs(c->fTracks[2])) max++;
2588   }
2589   if (max < Int_t(0.5*tail)) return -lab;
2590   
2591   return lab;
2592 }
2593
2594 //_____________________________________________________________________________
2595 void AliTPCtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const 
2596 {
2597   //-----------------------------------------------------------------
2598   // This function returns reconstructed track momentum in the global system.
2599   //
2600   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2601   //-----------------------------------------------------------------
2602   Double_t pt=TMath::Abs(GetPt()); // GeV/c
2603   Double_t r=x(2)*fX-x(3);
2604   Double_t y0=x(0) + sqrt(1.- r*r)/x(2);
2605   px=-pt*(x(0)-y0)*x(2);    //cos(phi);
2606   py=-pt*(x(3)-fX*x(2));   //sin(phi);
2607   pz=pt*x(4);
2608   Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha);
2609   py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha);
2610   px=tmp;  
2611 }
2612
2613 //_____________________________________________________________________________
2614 Double_t AliTPCtrack::GetdEdX(Double_t low, Double_t up) const {
2615   //-----------------------------------------------------------------
2616   // This funtion calculates dE/dX within the "low" and "up" cuts.
2617   //
2618   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2619   //-----------------------------------------------------------------
2620   Int_t ncl=fClusters.GetEntriesFast();
2621   Int_t n=0;
2622   Double_t *q=new Double_t[ncl];
2623   Int_t i;
2624   for (i=1; i<ncl; i++) { //Shall I think of this "i=1" ? (I.Belikov)
2625      AliTPCcluster *cl=(AliTPCcluster*)(fClusters.UncheckedAt(i));
2626      q[n++]=TMath::Abs(cl->fQ)/cl->fdEdX;
2627      if (cl->fSector<36) q[n-1]*=1.1;
2628   }
2629
2630   //stupid sorting
2631   Int_t swap;
2632   do {
2633     swap=0;
2634     for (i=0; i<n-1; i++) {
2635       if (q[i]<=q[i+1]) continue;
2636       Double_t tmp=q[i]; q[i]=q[i+1]; q[i+1]=tmp;
2637       swap++;
2638     }
2639   } while (swap);
2640
2641   Int_t nl=Int_t(low*n), nu=Int_t(up *n);
2642   Double_t dedx=0.;
2643   for (i=nl; i<=nu; i++) dedx += q[i];
2644   dedx /= (nu-nl+1);
2645   return dedx;
2646 }
2647
2648 //_________________________________________________________________________
2649 //
2650 // Classes for internal tracking use
2651 //_________________________________________________________________________
2652 void AliTPCRow::InsertCluster(const AliTPCcluster* c) {
2653   //-----------------------------------------------------------------------
2654   // Insert a cluster into this pad row in accordence with its y-coordinate
2655   //
2656   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2657   //-----------------------------------------------------------------------
2658   if (num_of_clusters==MAX_CLUSTER_PER_ROW) {
2659     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
2660   }
2661   if (num_of_clusters==0) {clusters[num_of_clusters++]=c; return;}
2662   Int_t i=Find(c->fY);
2663   memmove(clusters+i+1 ,clusters+i,(num_of_clusters-i)*sizeof(AliTPCcluster*));
2664   clusters[i]=c; num_of_clusters++;
2665 }
2666 //___________________________________________________________________
2667
2668 Int_t AliTPCRow::Find(Double_t y) const {
2669   //-----------------------------------------------------------------------
2670   // Return the index of the nearest cluster 
2671   //
2672   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
2673   //-----------------------------------------------------------------------
2674   if (y <= clusters[0]->fY) return 0;
2675   if (y > clusters[num_of_clusters-1]->fY) return num_of_clusters;
2676   Int_t b=0, e=num_of_clusters-1, m=(b+e)/2;
2677   for (; b<e; m=(b+e)/2) {
2678     if (y > clusters[m]->fY) b=m+1;
2679     else e=m; 
2680   }
2681   return m;
2682 }
2683 //________________________________________________________________________
2684