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