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