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="gif/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 "GParticle.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   AliMC* pMC = AliMC::GetMC();
 550
 551   Int_t ISXFLD=gAlice->Field()->Integ();
 552   Float_t SXMGMX=gAlice->Field()->Max();
 553
 554   Float_t absl, radl, a, d, z;
 555   Float_t dg;
 556   Float_t x0ne;
 557   Float_t buf[1];
 558   Int_t nbuf;
 559
 560   // --- Methane (CH4) ---
 561   Float_t am[2] = { 12.,1. };
 562   Float_t zm[2] = { 6.,1. };
 563   Float_t wm[2] = { 1.,4. };
 564   Float_t dm    = 7.17e-4;
 565   // --- The Neon CO2 90/10 mixture ---
 566   Float_t ag[2] = { 20.18 };
 567   Float_t zg[2] = { 10. };
 568   Float_t wg[2] = { .8,.2 };
 569   Float_t dne   = 9e-4;        // --- Neon density in g/cm3 ---
 570
 571   // --- Mylar (C5H4O2) ---
 572   Float_t amy[3] = { 12.,1.,16. };
 573   Float_t zmy[3] = { 6.,1.,8. };
 574   Float_t wmy[3] = { 5.,4.,2. };
 575   Float_t dmy    = 1.39;
 576   // --- CO2 ---
 577   Float_t ac[2] = { 12.,16. };
 578   Float_t zc[2] = { 6.,8. };
 579   Float_t wc[2] = { 1.,2. };
 580   Float_t dc    = .001977;
 581   // --- Carbon density and radiation length ---
 582   Float_t densc = 2.265;
 583   Float_t radlc = 18.8;
 584   // --- Silicon ---
 585   Float_t asi   = 28.09;
 586   Float_t zsi   = 14.;
 587   Float_t desi  = 2.33;
 588   Float_t radsi = 9.36;
 589
 590   // --- Define the various materials for GEANT ---
 591   AliMaterial(0, "Al $", 26.98, 13., 2.7, 8.9, 37.2);
 592   x0ne = 28.94 / dne;
 593   AliMaterial(1, "Ne $", 20.18, 10., dne, x0ne, 999.);
 594
 595   // --  Methane, defined by the proportions of atoms
 596
 597   AliMixture(2, "Methane$", am, zm, dm, -2, wm);
 598
 599   // --- CO2, defined by the proportion of atoms
 600
 601   AliMixture(7, "CO2$", ac, zc, dc, -2, wc);
 602
 603   // --  Get A,Z etc. for CO2
 604
 605   char namate[21];
 606   pMC->Gfmate((*fIdmate)[7], namate, a, z, d, radl, absl, buf, nbuf);
 607   ag[1] = a;
 608   zg[1] = z;
 609   dg = dne * .9 + dc * .1;
 610
 611   // --  Create Ne/CO2 90/10 mixture
 612
 613   AliMixture(3, "Gas-mixt $", ag, zg, dg, 2, wg);
 614   AliMixture(4, "Gas-mixt $", ag, zg, dg, 2, wg);
 615
 616   AliMaterial(5, "G10$", 20., 10., 1.7, 19.4, 999.);
 617   AliMixture(6, "Mylar$", amy, zmy, dmy, -3, wmy);
 618
 619   a = ac[0];
 620   z = zc[0];
 621   AliMaterial(8, "Carbon", a, z, densc, radlc, 999.);
 622
 623   AliMaterial(9, "Silicon", asi, zsi, desi, radsi, 999.);
 624   AliMaterial(99, "Air$", 14.61, 7.3, .001205, 30420., 67500.);
 625
 626   AliMedium(400, "Al wall$",  0, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1);
 627   AliMedium(402, "Gas mix1$", 3, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
 628   AliMedium(403, "Gas mix2$", 3, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
 629   AliMedium(404, "Gas mix3$", 4, 1, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
 630   AliMedium(405, "G10 pln$",  5, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1 );
 631   AliMedium(406, "Mylar  $",  6, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
 632   AliMedium(407, "CO2    $",  7, 0, ISXFLD, SXMGMX, 10., .01,.1, .01,  .01);
 633   AliMedium(408, "Carbon $",  8, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1 );
 634   AliMedium(409, "Silicon$",  9, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1 );
 635   AliMedium(499, "Air gap$", 99, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1 );
 636 }
 637
 638 //_____________________________________________________________________________
 639 struct Bin {
 640    const AliTPCdigit *dig;
 641    int idx;
 642    Bin() {dig=0; idx=-1;}
 643 };
 644
 645 struct PreCluster : public AliTPCcluster {
 646    const AliTPCdigit* summit;
 647    int idx;
 648    int cut;
 649    int npeaks;
 650    PreCluster();
 651 };
 652 PreCluster::PreCluster() : AliTPCcluster() {cut=npeaks=0;}
 653
 654
 655 //_____________________________________________________________________________
 656 static void FindCluster(int i, int j, Bin bins[][MAXTPCTBK+2], PreCluster &c)
 657 {
 658   //
 659   // Find clusters
 660   //
 661   Bin& b=bins[i][j];
 662   double q=double(b.dig->fSignal);
 663
 664   if (q<0) { // digit is at the edge of the pad row
 665     q=-q;
 666     c.cut=1;
 667   }
 668   if (b.idx >= 0 && b.idx != c.idx) {
 669     c.idx=b.idx;
 670     c.npeaks++;
 671   }
 672
 673   if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
 674
 675   c.fY += i*q;
 676   c.fZ += j*q;
 677   c.fSigmaY2 += i*i*q;
 678   c.fSigmaZ2 += j*j*q;
 679   c.fQ += q;
 680
 681   b.dig = 0;  b.idx = c.idx;
 682
 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   if (bins[i+1][j].dig) FindCluster(i+1,j,bins,c);
 686   if (bins[i][j+1].dig) FindCluster(i,j+1,bins,c);
 687
 688 }
 689
 690 //_____________________________________________________________________________
 691 void AliTPC::Digits2Clusters()
 692 {
 693   //
 694   // simple TPC cluster finder from digits.
 695   //
 696   //
 697   AliTPCParam * fTPCParam = &(fDigParam->GetParam());
 698
 699   const Int_t MAX_PAD=200+2, MAX_BUCKET=MAXTPCTBK+2;
 700   const Int_t Q_min=60;
 701   const Int_t THRESHOLD=20;
 702
 703   TTree *t=(TTree*)gDirectory->Get("TreeD0_Param1");
 704   t->GetBranch("Digits")->SetAddress(&fDigits);
 705   Int_t sectors_by_rows=(Int_t)t->GetEntries();
 706
 707   int ncls=0;
 708
 709   for (Int_t n=0; n<sectors_by_rows; n++) {
 710     if (!t->GetEvent(n)) continue;
 711     Bin bins[MAX_PAD][MAX_BUCKET];
 712     AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0);
 713     Int_t nsec=dig->fSector, nrow=dig->fPadRow;
 714     Int_t ndigits=fDigits->GetEntriesFast();
 715
 716     int npads;  int sign_z;
 717     if (nsec<25) {
 718       sign_z=(nsec<13) ? 1 : -1;
 719       npads=fTPCParam->GetNPadsLow(nrow-1);
 720     } else {
 721       sign_z=(nsec<49) ? 1 : -1;
 722       npads=fTPCParam->GetNPadsUp(nrow-1);
 723     }
 724
 725     int ndig;
 726     for (ndig=0; ndig<ndigits; ndig++) {
 727       dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
 728       int i=dig->fPad, j=dig->fTime;
 729       if (dig->fSignal >= THRESHOLD) bins[i][j].dig=dig;
 730       if (i==1 || i==npads || j==1 || j==MAXTPCTBK) dig->fSignal*=-1;
 731     }
 732
 733     int ncl=0;
 734     int i,j;
 735
 736     for (i=1; i<MAX_PAD-1; i++) {
 737       for (j=1; j<MAX_BUCKET-1; j++) {
 738         if (bins[i][j].dig == 0) continue;
 739         PreCluster c; c.summit=bins[i][j].dig; c.idx=ncls;
 740         FindCluster(i, j, bins, c);
 741         c.fY /= c.fQ;
 742         c.fZ /= c.fQ;
 743
 744         double s2 = c.fSigmaY2/c.fQ - c.fY*c.fY;
 745         c.fSigmaY2 = s2 + 1./12.;
 746         c.fSigmaY2 *= fTPCParam->GetPadPitchWidth()*
 747                       fTPCParam->GetPadPitchWidth();
 748         if (s2 != 0.) c.fSigmaY2 *= 0.022*8*4;
 749
 750         s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ;
 751         c.fSigmaZ2 = s2 + 1./12.;
 752         c.fSigmaZ2 *= fTPCParam->GetZWidth()*fTPCParam->GetZWidth();
 753         if (s2 != 0.) c.fSigmaZ2 *= 0.068*4*4;
 754
 755         c.fY = (c.fY - 0.5 - 0.5*npads)*fTPCParam->GetPadPitchWidth();
 756         c.fZ = fTPCParam->GetZWidth()*(c.fZ+1);
 757         c.fZ -= 3.*fTPCParam->GetZSigma(); // PASA delay
 758         c.fZ = sign_z*(z_end - c.fZ);
 759         //c.fZ += 0.023;
 760         c.fSector=nsec;
 761         c.fPadRow=nrow;
 762         c.fTracks[0]=c.summit->fTracks[0];
 763         c.fTracks[1]=c.summit->fTracks[1];
 764         c.fTracks[2]=c.summit->fTracks[2];
 765
 766         if (c.cut) {
 767           c.fSigmaY2 *= 25.;
 768           c.fSigmaZ2 *= 4.;
 769         }
 770
 771         AddCluster(c); ncls++; ncl++;
 772       }
 773     }
 774
 775     for (ndig=0; ndig<ndigits; ndig++) {
 776       dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
 777       if (TMath::Abs(dig->fSignal) >= 0)
 778         bins[dig->fPad][dig->fTime].dig=dig;
 779     }
 780
 781     for (i=1; i<MAX_PAD-1; i++) {
 782       for (j=1; j<MAX_BUCKET-1; j++) {
 783         if (bins[i][j].dig == 0) continue;
 784         PreCluster c; c.summit=bins[i][j].dig; c.idx=ncls;
 785         FindCluster(i, j, bins, c);
 786         if (c.fQ <= Q_min) continue; //noise cluster
 787         if (c.npeaks>1) continue;    //overlapped cluster
 788         c.fY /= c.fQ;
 789         c.fZ /= c.fQ;
 790
 791         double s2 = c.fSigmaY2/c.fQ - c.fY*c.fY;
 792         c.fSigmaY2 = s2 + 1./12.;
 793         c.fSigmaY2 *= fTPCParam->GetPadPitchWidth()*
 794                       fTPCParam->GetPadPitchWidth();
 795         if (s2 != 0.) c.fSigmaY2 *= 0.022*4*0.6*4;
 796
 797         s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ;
 798         c.fSigmaZ2 = s2 + 1./12.;
 799         c.fSigmaZ2 *= fTPCParam->GetZWidth()*fTPCParam->GetZWidth();
 800         if (s2 != 0.) c.fSigmaZ2 *= 0.068*4*0.4;
 801
 802         c.fY = (c.fY - 0.5 - 0.5*npads)*fTPCParam->GetPadPitchWidth();
 803         c.fZ = fTPCParam->GetZWidth()*(c.fZ+1);
 804         c.fZ -= 3.*fTPCParam->GetZSigma(); // PASA delay
 805         c.fZ = sign_z*(z_end - c.fZ);
 806         //c.fZ += 0.023;
 807         c.fSector=nsec;
 808         c.fPadRow=nrow;
 809         c.fTracks[0]=c.summit->fTracks[0];
 810         c.fTracks[1]=c.summit->fTracks[1];
 811         c.fTracks[2]=c.summit->fTracks[2];
 812
 813         if (c.cut) {
 814           c.fSigmaY2 *= 25.;
 815           c.fSigmaZ2 *= 4.;
 816         }
 817
 818         if (c.npeaks==0) {AddCluster(c); ncls++; ncl++;}
 819         else {
 820           new ((*fClusters)[c.idx]) AliTPCcluster(c);
 821         }
 822       }
 823     }
 824
 825     cerr<<"sector, row, digits, clusters: "
 826         <<nsec<<' '<<nrow<<' '<<ndigits<<' '<<ncl<<"                  \r";
 827
 828     fDigits->Clear();
 829
 830   }
 831 }
 832
 833 //_____________________________________________________________________________
 834 void AliTPC::ElDiff(Float_t *xyz)
 835 {
 836   //--------------------------------------------------
 837   // calculates the diffusion of a single electron
 838   //--------------------------------------------------
 839
 840   //-----------------------------------------------------------------
 841   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
 842   //-----------------------------------------------------------------
 843   AliTPCParam * fTPCParam = &(fDigParam->GetParam());
 844   Float_t driftl;
 845
 846   Float_t z0=xyz[2];
 847
 848   driftl=z_end-TMath::Abs(xyz[2]);
 849
 850   if(driftl<0.01) driftl=0.01;
 851
 852   // check the attachment
 853
 854   driftl=TMath::Sqrt(driftl);
 855
 856   //  Float_t sig_t = driftl*diff_t;
 857   //Float_t sig_l = driftl*diff_l;
 858   Float_t sig_t = driftl*fTPCParam->GetDiffT();
 859   Float_t sig_l = driftl*fTPCParam->GetDiffL();
 860   xyz[0]=gRandom->Gaus(xyz[0],sig_t);
 861   xyz[1]=gRandom->Gaus(xyz[1],sig_t);
 862   xyz[2]=gRandom->Gaus(xyz[2],sig_l);
 863
 864   if (TMath::Abs(xyz[2])>z_end){
 865     xyz[2]=TMath::Sign(z_end,z0);
 866   }
 867   if(xyz[2]*z0 < 0.){
 868     Float_t eps = 0.0001;
 869     xyz[2]=TMath::Sign(eps,z0);
 870   }
 871 }
 872
 873 //_____________________________________________________________________________
 874 void AliTPC::Hits2Clusters()
 875 {
 876   //--------------------------------------------------------
 877   // TPC simple cluster generator from hits
 878   // obtained from the TPC Fast Simulator
 879   // The point errors are taken from the parametrization
 880   //--------------------------------------------------------
 881
 882   //-----------------------------------------------------------------
 883   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
 884   //-----------------------------------------------------------------
 885   AliTPCParam * fTPCParam = &(fDigParam->GetParam());
 886   Float_t sigma_rphi,sigma_z,cl_rphi,cl_z;
 887   //
 888   GParticle *particle; // pointer to a given particle
 889   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 890   TClonesArray *Particles; //pointer to the particle list
 891   Int_t sector,nhits;
 892   Int_t ipart;
 893   Float_t xyz[5];
 894   Float_t pl,pt,tanth,rpad,ratio;
 895   Float_t cph,sph;
 896
 897   //---------------------------------------------------------------
 898   //  Get the access to the tracks
 899   //---------------------------------------------------------------
 900
 901   TTree *TH = gAlice->TreeH();
 902   Stat_t ntracks = TH->GetEntries();
 903
 904   //------------------------------------------------------------
 905   // Loop over all sectors (72 sectors)
 906   // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0
 907   // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0
 908   //
 909   // First cluster for sector 1 starts at "0"
 910   //------------------------------------------------------------
 911
 912
 913   fClustersIndex[0] = 0;
 914
 915   //
 916   for(Int_t isec=1;isec<fNsectors+1;isec++){
 917     //MI change
 918     fTPCParam->AdjustAngles(isec,cph,sph);
 919
 920     //------------------------------------------------------------
 921     // Loop over tracks
 922     //------------------------------------------------------------
 923
 924     for(Int_t track=0;track<ntracks;track++){
 925       ResetHits();
 926       TH->GetEvent(track);
 927       //
 928       //  Get number of the TPC hits and a pointer
 929       //  to the particles
 930       //
 931       nhits=fHits->GetEntriesFast();
 932       Particles=gAlice->Particles();
 933       //
 934       // Loop over hits
 935       //
 936       for(Int_t hit=0;hit<nhits;hit++){
 937         tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
 938         sector=tpcHit->fSector; // sector number
 939         if(sector != isec) continue; //terminate iteration
 940         ipart=tpcHit->fTrack;
 941         particle=(GParticle*)Particles->UncheckedAt(ipart);
 942         pl=particle->GetPz();
 943         pt=particle->GetPT();
 944         if(pt < 1.e-9) pt=1.e-9;
 945         tanth=pl/pt;
 946         tanth = TMath::Abs(tanth);
 947         rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY);
 948         ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
 949
 950         //   space-point resolutions
 951
 952         sigma_rphi=SigmaY2(rpad,tanth,pt);
 953         sigma_z   =SigmaZ2(rpad,tanth   );
 954
 955         //   cluster widths
 956
 957         cl_rphi=ac_rphi-bc_rphi*rpad*tanth+cc_rphi*ratio*ratio;
 958         cl_z=ac_z-bc_z*rpad*tanth+cc_z*tanth*tanth;
 959
 960         // temporary protection
 961
 962         if(sigma_rphi < 0.) sigma_rphi=0.4e-3;
 963         if(sigma_z < 0.) sigma_z=0.4e-3;
 964         if(cl_rphi < 0.) cl_rphi=2.5e-3;
 965         if(cl_z < 0.) cl_z=2.5e-5;
 966
 967         //
 968
 969         //
 970         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
 971         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
 972         //
 973         Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph;
 974         Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph;
 975         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigma_rphi));   // y
 976         Double_t alpha=(sector<25) ? alpha_low : alpha_up;
 977         if (TMath::Abs(xyz[0]/xprim) > TMath::Tan(0.5*alpha)) xyz[0]=yprim;
 978         xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigma_z)); // z
 979         if (TMath::Abs(xyz[1]) > 250) xyz[1]=tpcHit->fZ;
 980         xyz[2]=tpcHit->fQ;                                     // q
 981         xyz[3]=sigma_rphi;                                     // fSigmaY2
 982         xyz[4]=sigma_z;                                        // fSigmaZ2
 983
 984         //find row number
 985         //MI we must change
 986         Int_t row = fTPCParam->GetPadRow(sector,xprim) ;
 987         // and finally add the cluster
 988         Int_t tracks[5]={tpcHit->fTrack, -1, -1, sector, row+1};
 989         AddCluster(xyz,tracks);
 990
 991       } // end of loop over hits
 992     }   // end of loop over tracks
 993
 994     fClustersIndex[isec] = fNclusters; // update clusters index
 995
 996   } // end of loop over sectors
 997
 998   fClustersIndex[fNsectors]--; // set end of the clusters buffer
 999
1000 }
1001
1002
1003 void AliTPC::Hits2Digits()
1004 {
1005
1006  //----------------------------------------------------
1007   // Loop over all sectors (72 sectors)
1008   // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0
1009   // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0
1010   //----
1011   for(Int_t isec=1;isec<fNsectors+1;isec++)  Hits2DigitsSector(isec);
1012 }
1013
1014
1015 //_____________________________________________________________________________
1016 void AliTPC::Hits2DigitsSector(Int_t isec)
1017 {
1018   //-------------------------------------------------------------------
1019   // TPC conversion from hits to digits.
1020   //-------------------------------------------------------------------
1021
1022   //-----------------------------------------------------------------
1023   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1024   //-----------------------------------------------------------------
1025
1026   //-------------------------------------------------------
1027   //  Get the access to the track hits
1028   //-------------------------------------------------------
1029
1030   AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1031   TTree *TH = gAlice->TreeH(); // pointer to the hits tree
1032   Stat_t ntracks = TH->GetEntries();
1033
1034   if( ntracks > 0){
1035
1036   //-------------------------------------------
1037   //  Only if there are any tracks...
1038   //-------------------------------------------
1039
1040
1041     // TObjArrays for three neighouring pad-rows
1042
1043     TObjArray **rowTriplet = new TObjArray* [3];
1044
1045     // TObjArray-s for each pad-row
1046
1047     TObjArray **row;
1048
1049     for (Int_t trip=0;trip<3;trip++){
1050       rowTriplet[trip]=new TObjArray;
1051     }
1052
1053
1054
1055       printf("*** Processing sector number %d ***\n",isec);
1056
1057       Int_t nrows =fTPCParam->GetNRow(isec);
1058
1059       row= new TObjArray* [nrows];
1060
1061       MakeSector(isec,nrows,TH,ntracks,row);
1062
1063       //--------------------------------------------------------
1064       //   Digitize this sector, row by row
1065       //   row[i] is the pointer to the TObjArray of TVectors,
1066       //   each one containing electrons accepted on this
1067       //   row, assigned into tracks
1068       //--------------------------------------------------------
1069
1070       Int_t i;
1071
1072       for (i=0;i<nrows;i++){
1073
1074         // Triplets for i = 0 and i=1 are identical!
1075         // The same for i = nrows-1 and nrows!
1076
1077         if(i != 1 && i != nrows-1){
1078            MakeTriplet(i,rowTriplet,row);
1079          }
1080
1081             DigitizeRow(i,isec,rowTriplet);
1082
1083             fDigParam->Fill();
1084
1085             Int_t ndig=fDigParam->GetArray()->GetEntriesFast();
1086
1087             printf("*** Sector, row, digits %d %d %d ***\n",isec,i,ndig);
1088
1089             ResetDigits(); // reset digits for this row after storing them
1090
1091        } // end of the sector digitization
1092
1093        // delete the last triplet
1094
1095        for (i=0;i<3;i++) rowTriplet[i]->Delete();
1096
1097        delete [] row; // delete the array of pointers to TObjArray-s
1098
1099   } // ntracks >0
1100 } // end of Hits2Digits
1101 //_____________________________________________________________________________
1102 void AliTPC::MakeTriplet(Int_t row,
1103                          TObjArray **rowTriplet, TObjArray **prow)
1104 {
1105   //------------------------------------------------------------------
1106   //  Makes the "triplet" of the neighbouring pad-row for the
1107   //  digitization including the cross-talk between the pad-rows
1108   //------------------------------------------------------------------
1109
1110   //-----------------------------------------------------------------
1111   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1112   //-----------------------------------------------------------------
1113
1114   AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1115   Float_t gasgain = fTPCParam->GetGasGain();
1116   Int_t nTracks[3];
1117
1118   Int_t nElements,nElectrons;
1119
1120   TVector *pv;
1121   TVector *tv;
1122
1123   //-------------------------------------------------------------------
1124   // pv is an "old" track, i.e. label + triplets of (x,y,z)
1125   // for each electron
1126   //
1127   //-------------------------------------------------------------------
1128
1129
1130   Int_t i1,i2;
1131   Int_t nel,nt;
1132
1133   if(row == 0 || row == 1){
1134
1135   // create entire triplet for the first AND the second row
1136
1137     nTracks[0] = prow[0]->GetEntries();
1138     nTracks[1] = prow[1]->GetEntries();
1139     nTracks[2] = prow[2]->GetEntries();
1140
1141     for(i1=0;i1<3;i1++){
1142       nt = nTracks[i1]; // number of tracks for this row
1143
1144       for(i2=0;i2<nt;i2++){
1145         pv = (TVector*)prow[i1]->At(i2);
1146         TVector &v1 = *pv;
1147         nElements = pv->GetNrows();
1148         nElectrons = (nElements-1)/3;
1149
1150         tv = new TVector(4*nElectrons+1); // create TVector for a modified track
1151         TVector &v2 = *tv;
1152         v2(0)=v1(0); //track label
1153
1154         for(nel=0;nel<nElectrons;nel++){
1155           Int_t idx1 = nel*3;
1156           Int_t idx2 = nel*4;
1157           // Avalanche, including fluctuations
1158           Int_t aval = (Int_t) (-gasgain*TMath::Log(gRandom->Rndm()));
1159           v2(idx2+1)= v1(idx1+1);
1160           v2(idx2+2)= v1(idx1+2);
1161           v2(idx2+3)= v1(idx1+3);
1162           v2(idx2+4)= (Float_t)aval; // in number of electrons!
1163         } // end of loop over electrons
1164         //
1165         //  Add this track to a row
1166         //
1167
1168         rowTriplet[i1]->Add(tv);
1169
1170
1171       } // end of loop over tracks for this row
1172
1173       prow[i1]->Delete(); // remove "old" tracks
1174       delete prow[i1]; // delete  TObjArray for this row
1175       prow[i1]=0; // set pointer to NULL
1176
1177     } // end of loop over row triplets
1178
1179
1180   }
1181   else{
1182
1183     rowTriplet[0]->Delete(); // remove old lower row
1184
1185     nTracks[0]=rowTriplet[1]->GetEntries(); // previous middle row
1186     nTracks[1]=rowTriplet[2]->GetEntries(); // previous upper row
1187     nTracks[2]=prow[row+1]->GetEntries(); // next row
1188
1189
1190     //-------------------------------------------
1191     //  shift new tracks downwards
1192     //-------------------------------------------
1193
1194     for(i1=0;i1<nTracks[0];i1++){
1195       pv=(TVector*)rowTriplet[1]->At(i1);
1196       rowTriplet[0]->Add(pv);
1197     }
1198
1199     rowTriplet[1]->Clear(); // leave tracks on the heap!!!
1200
1201     for(i1=0;i1<nTracks[1];i1++){
1202       pv=(TVector*)rowTriplet[2]->At(i1);
1203       rowTriplet[1]->Add(pv);
1204     }
1205
1206     rowTriplet[2]->Clear(); // leave tracks on the heap!!!
1207
1208     //---------------------------------------------
1209     //  Create new upper row
1210     //---------------------------------------------
1211
1212
1213
1214     for(i1=0;i1<nTracks[2];i1++){
1215         pv = (TVector*)prow[row+1]->At(i1);
1216         TVector &v1 = *pv;
1217         nElements = pv->GetNrows();
1218         nElectrons = (nElements-1)/3;
1219
1220         tv = new TVector(4*nElectrons+1); // create TVector for a modified track
1221         TVector &v2 = *tv;
1222         v2(0)=v1(0); //track label
1223
1224         for(nel=0;nel<nElectrons;nel++){
1225
1226           Int_t idx1 = nel*3;
1227           Int_t idx2 = nel*4;
1228           // Avalanche, including fluctuations
1229           Int_t aval = (Int_t)
1230             (-gasgain*TMath::Log(gRandom->Rndm()));
1231
1232           v2(idx2+1)= v1(idx1+1);
1233           v2(idx2+2)= v1(idx1+2);
1234           v2(idx2+3)= v1(idx1+3);
1235           v2(idx2+4)= (Float_t)aval; // in number of electrons!
1236         } // end of loop over electrons
1237
1238           rowTriplet[2]->Add(tv);
1239
1240     } // end of loop over tracks
1241
1242     prow[row+1]->Delete(); // delete tracks for this row
1243     delete prow[row+1];  // delete TObjArray for this row
1244     prow[row+1]=0; // set a pointer to NULL
1245
1246   }
1247
1248 }  // end of MakeTriplet
1249 //_____________________________________________________________________________
1250 void AliTPC::ExB(Float_t *xyz)
1251 {
1252   //------------------------------------------------
1253   //  ExB at the wires and wire number calulation
1254   //------------------------------------------------
1255
1256   //-----------------------------------------------------------------
1257   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1258   //-----------------------------------------------------------------
1259    AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1260
1261    Float_t x1=xyz[0];
1262    fTPCParam->GetWire(x1);        //calculate nearest wire position
1263    Float_t dx=xyz[0]-x1;
1264    xyz[1]+=dx*fTPCParam->GetOmegaTau();
1265
1266 } // end of ExB
1267 //_____________________________________________________________________________
1268 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rowTriplet)
1269 {
1270   //-----------------------------------------------------------
1271   // Single row digitization, coupling from the neighbouring
1272   // rows taken into account
1273   //-----------------------------------------------------------
1274
1275   //-----------------------------------------------------------------
1276   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1277   //-----------------------------------------------------------------
1278
1279
1280   AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1281   Float_t chipgain= fTPCParam->GetChipGain();
1282   Float_t zerosup = fTPCParam->GetZeroSup();
1283   Int_t nrows =fTPCParam->GetNRow(isec);
1284
1285   Int_t nTracks[3];
1286   Int_t n_of_pads[3];
1287   Int_t IndexRange[4];
1288   Int_t i1;
1289   Int_t iFlag;
1290
1291   //
1292   // iFlag = 0 -> inner row, iFlag = 1 -> the middle one, iFlag = 2 -> the outer one
1293   //
1294
1295   nTracks[0]=rowTriplet[0]->GetEntries(); // lower row
1296   nTracks[1]=rowTriplet[1]->GetEntries(); // middle row
1297   nTracks[2]=rowTriplet[2]->GetEntries(); // upper row
1298
1299
1300   if(irow == 0){
1301     iFlag=0;
1302     n_of_pads[0]=fTPCParam->GetNPads(isec,0);
1303     n_of_pads[1]=fTPCParam->GetNPads(isec,1);
1304   }
1305   else if(irow == nrows-1){
1306      iFlag=2;
1307      n_of_pads[1]=fTPCParam->GetNPads(isec,irow-1);
1308      n_of_pads[2]=fTPCParam->GetNPads(isec,irow);
1309   }
1310   else {
1311     iFlag=1;
1312     for(i1=0;i1<3;i1++){
1313        n_of_pads[i1]=fTPCParam->GetNPads(isec,irow-1+i1);
1314     }
1315   }
1316
1317   //
1318   //  Integrated signal for this row
1319   //  and a single track signal
1320   //
1321
1322   TMatrix *m1   = new TMatrix(1,n_of_pads[iFlag],1,MAXTPCTBK); // integrated
1323   TMatrix *m2   = new TMatrix(1,n_of_pads[iFlag],1,MAXTPCTBK); // single
1324
1325   //
1326
1327   TMatrix &Total  = *m1;
1328
1329   //  Array of pointers to the label-signal list
1330
1331   Int_t NofDigits = n_of_pads[iFlag]*MAXTPCTBK; // number of digits for this row
1332
1333   Float_t  **pList = new Float_t* [NofDigits];
1334
1335   Int_t lp;
1336
1337   for(lp=0;lp<NofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1338
1339   //
1340   //  Straight signal and cross-talk, cross-talk is integrated over all
1341   //  tracks and added to the signal at the very end
1342   //
1343
1344
1345   for (i1=0;i1<nTracks[iFlag];i1++){
1346
1347     m2->Zero(); // clear single track signal matrix
1348
1349     Float_t TrackLabel =
1350       GetSignal(rowTriplet[iFlag],i1,n_of_pads[iFlag],m2,m1,IndexRange);
1351
1352     GetList(TrackLabel,n_of_pads[iFlag],m2,IndexRange,pList);
1353
1354   }
1355
1356   //
1357   //  Cross talk from the neighbouring pad-rows
1358   //
1359
1360   TMatrix *m3 =  new TMatrix(1,n_of_pads[iFlag],1,MAXTPCTBK); // cross-talk
1361
1362   TMatrix &Cross = *m3;
1363
1364   if(iFlag == 0){
1365
1366     // cross-talk from the outer row only (first pad row)
1367
1368     GetCrossTalk(0,rowTriplet[1],nTracks[1],n_of_pads,m3);
1369
1370   }
1371   else if(iFlag == 2){
1372
1373     // cross-talk from the inner row only (last pad row)
1374
1375     GetCrossTalk(2,rowTriplet[1],nTracks[1],n_of_pads,m3);
1376
1377   }
1378   else{
1379
1380     // cross-talk from both inner and outer rows
1381
1382     GetCrossTalk(3,rowTriplet[0],nTracks[0],n_of_pads,m3); // inner
1383     GetCrossTalk(4,rowTriplet[2],nTracks[2],n_of_pads,m3); //outer
1384   }
1385
1386   Total += Cross; // add the cross-talk
1387
1388   //
1389   //  Convert analog signal to ADC counts
1390   //
1391
1392   Int_t tracks[3];
1393   Int_t digits[5];
1394
1395
1396   for(Int_t ip=1;ip<n_of_pads[iFlag]+1;ip++){
1397     for(Int_t it=1;it<MAXTPCTBK+1;it++){
1398
1399       Float_t q = Total(ip,it);
1400
1401       Int_t gi =(it-1)*n_of_pads[iFlag]+ip-1; // global index
1402
1403       q = gRandom->Gaus(q,fTPCParam->GetNoise()); // apply noise
1404       q *= (q_el*1.e15); // convert to fC
1405       q *= chipgain; // convert to mV
1406       q *= (adc_sat/dyn_range); // convert to ADC counts
1407
1408       if(q <zerosup) continue; // do not fill zeros
1409       if(q > adc_sat) q = adc_sat;  // saturation
1410
1411       //
1412       //  "real" signal or electronic noise (list = -1)?
1413       //
1414
1415       for(Int_t j1=0;j1<3;j1++){
1416         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1;
1417       }
1418
1419       digits[0]=isec;
1420       digits[1]=irow+1;
1421       digits[2]=ip;
1422       digits[3]=it;
1423       digits[4]= (Int_t)q;
1424
1425       //  Add this digit
1426
1427       AddDigit(tracks,digits);
1428
1429     } // end of loop over time buckets
1430   }  // end of lop over pads
1431
1432   //
1433   //  This row has been digitized, delete nonused stuff
1434   //
1435
1436   for(lp=0;lp<NofDigits;lp++){
1437     if(pList[lp]) delete [] pList[lp];
1438   }
1439
1440   delete [] pList;
1441
1442   delete m1;
1443   delete m2;
1444   delete m3;
1445
1446 } // end of DigitizeRow
1447 //_____________________________________________________________________________
1448 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, Int_t np, TMatrix *m1, TMatrix *m2,
1449                           Int_t *IndexRange)
1450 {
1451
1452   //---------------------------------------------------------------
1453   //  Calculates 2-D signal (pad,time) for a single track,
1454   //  returns a pointer to the signal matrix and the track label
1455   //  No digitization is performed at this level!!!
1456   //---------------------------------------------------------------
1457
1458   //-----------------------------------------------------------------
1459   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1460   //-----------------------------------------------------------------
1461
1462   TVector *tv;
1463   AliTPCParam * fTPCParam = &(fDigParam->GetParam());
1464   AliTPCPRF2D * fPRF2D = &(fDigParam->GetPRF2D());
1465   AliTPCRF1D  * fRF    = &(fDigParam->GetRF());
1466
1467   //to make the code faster we put parameters  to the stack
1468
1469   Float_t zwidth  = fTPCParam->GetZWidth();
1470   Float_t zwidthm1  =1./fTPCParam->GetZWidth();
1471
1472   tv = (TVector*)p1->At(ntr); // pointer to a track
1473   TVector &v = *tv;
1474
1475   Float_t label = v(0);
1476
1477   Int_t CentralPad = (np+1)/2;
1478   Int_t PadNumber;
1479   Int_t nElectrons = (tv->GetNrows()-1)/4;
1480   Float_t range=((np-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); // pad range
1481   range -= 0.5; // dead zone, 5mm from the edge, according to H.G. Fischer
1482
1483   Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
1484
1485
1486   Float_t PadSignal[7]; // signal from a single electron
1487
1488   TMatrix &signal = *m1;
1489   TMatrix &total = *m2;
1490
1491
1492   IndexRange[0]=9999; // min pad
1493   IndexRange[1]=-1; // max pad
1494   IndexRange[2]=9999; //min time
1495   IndexRange[3]=-1; // max time
1496
1497   //
1498   //  Loop over all electrons
1499   //
1500
1501   for(Int_t nel=0; nel<nElectrons; nel++){
1502    Int_t idx=nel*4;
1503    Float_t xwire = v(idx+1);
1504    Float_t y = v(idx+2);
1505    Float_t z = v(idx+3);
1506
1507
1508    Float_t absy=TMath::Abs(y);
1509
1510    if(absy < 0.5*fTPCParam->GetPadPitchWidth()){
1511      PadNumber=CentralPad;
1512    }
1513    else if (absy < range){
1514      PadNumber=(Int_t) ((absy-0.5*fTPCParam->GetPadPitchWidth())/fTPCParam->GetPadPitchWidth() +1.);
1515      PadNumber=(Int_t) (TMath::Sign((Float_t)PadNumber, y)+CentralPad);
1516    }
1517    else continue; // electron out of pad-range , lost at the sector edge
1518
1519    Float_t aval = (absy<range-0.5) ? v(idx+4):v(idx+4)*IneffFactor;
1520
1521
1522    Float_t dist = y - (Float_t)(PadNumber-CentralPad)*fTPCParam->GetPadPitchWidth();
1523    for (Int_t i=0;i<7;i++)
1524      PadSignal[i]=fPRF2D->GetPRF(dist+(i-3)*fTPCParam->GetPadPitchWidth(),xwire)*aval;
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      // real pad range
1923
1924      Int_t  LeftPad = TMath::Max(1,PadNumber-3);
1925      Int_t  RightPad = TMath::Min(nPadsSignal,PadNumber+3);
1926
1927      Int_t pmin=LeftPad-PadNumber+3; // lower index of the pad_signal vector
1928      Int_t pmax=RightPad-PadNumber+3; // upper index
1929
1930
1931      Float_t z_drift = (z_end-z)*zwidthm1;
1932      Float_t z_offset = z_drift-(Int_t)z_drift;
1933      //distance to the centre of nearest time bin (in time bin units)
1934      Int_t FirstBucket = (Int_t)z_drift+1;
1935      // MI check it --time offset
1936      for (Int_t i2=0;i2<4;i2++){
1937        Int_t TrueTime = FirstBucket+i2; // current time bucket
1938        Float_t dz   = (Float_t(i2)+z_offset)*zwidth;
1939        Float_t ampl = fRF->GetRF(dz);
1940        if((TrueTime>MAXTPCTBK)) break; // beyond the time range
1941
1942
1943        // loop over pads, from pmin to pmax
1944
1945        for(Int_t i3=pmin;i3<pmax+1;i3++){
1946          Int_t TruePad = LeftPad+i3-pmin;
1947
1948          if(TruePad<nPadsDiff+1 || TruePad > nPadsSignal-nPadsDiff) continue;
1949
1950          TruePad -= nPadsDiff;
1951          signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!
1952
1953        } // end of loop over pads
1954      } // end of loop over time bins
1955
1956    } // end of loop over electrons
1957
1958  } // end of loop over tracks
1959
1960 } // end of GetCrossTalk
1961
1962
1963
1964 //_____________________________________________________________________________
1965 void AliTPC::Init()
1966 {
1967   //
1968   // Initialise TPC detector after definition of geometry
1969   //
1970   Int_t i;
1971   //
1972   printf("\n");
1973   for(i=0;i<35;i++) printf("*");
1974   printf(" TPC_INIT ");
1975   for(i=0;i<35;i++) printf("*");
1976   printf("\n");
1977   //
1978   for(i=0;i<80;i++) printf("*");
1979   printf("\n");
1980 }
1981
1982 //_____________________________________________________________________________
1983 void AliTPC::MakeBranch(Option_t* option)
1984 {
1985   //
1986   // Create Tree branches for the TPC.
1987   //
1988   Int_t buffersize = 4000;
1989   char branchname[10];
1990   sprintf(branchname,"%s",GetName());
1991
1992   AliDetector::MakeBranch(option);
1993
1994   char *D = strstr(option,"D");
1995
1996   if (fDigits   && gAlice->TreeD() && D) {
1997     gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
1998     printf("Making Branch %s for digits\n",branchname);
1999   }
2000
2001   char *R = strstr(option,"R");
2002
2003   if (fClusters && gAlice->TreeR() && R) {
2004     gAlice->TreeR()->Branch(branchname,&fClusters, buffersize);
2005     printf("Making Branch %s for Clusters\n",branchname);
2006   }
2007 }
2008
2009 //_____________________________________________________________________________
2010 void AliTPC::ResetDigits()
2011 {
2012   //
2013   // Reset number of digits and the digits array for this detector
2014   // reset clusters
2015   //
2016   fNdigits   = 0;
2017   //  if (fDigits)   fDigits->Clear();
2018   if (fDigParam->GetArray()!=0)  fDigParam->GetArray()->Clear();
2019   fNclusters = 0;
2020   if (fClusters) fClusters->Clear();
2021 }
2022
2023 //_____________________________________________________________________________
2024 void AliTPC::SetSecAL(Int_t sec)
2025 {
2026   //---------------------------------------------------
2027   // Activate/deactivate selection for lower sectors
2028   //---------------------------------------------------
2029
2030   //-----------------------------------------------------------------
2031   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2032   //-----------------------------------------------------------------
2033
2034   fSecAL = sec;
2035 }
2036
2037 //_____________________________________________________________________________
2038 void AliTPC::SetSecAU(Int_t sec)
2039 {
2040   //----------------------------------------------------
2041   // Activate/deactivate selection for upper sectors
2042   //---------------------------------------------------
2043
2044   //-----------------------------------------------------------------
2045   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2046   //-----------------------------------------------------------------
2047
2048   fSecAU = sec;
2049 }
2050
2051 //_____________________________________________________________________________
2052 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
2053 {
2054   //----------------------------------------
2055   // Select active lower sectors
2056   //----------------------------------------
2057
2058   //-----------------------------------------------------------------
2059   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2060   //-----------------------------------------------------------------
2061
2062   fSecLows[0] = s1;
2063   fSecLows[1] = s2;
2064   fSecLows[2] = s3;
2065   fSecLows[3] = s4;
2066   fSecLows[4] = s5;
2067   fSecLows[5] = s6;
2068 }
2069
2070 //_____________________________________________________________________________
2071 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
2072                        Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
2073                        Int_t s11 , Int_t s12)
2074 {
2075   //--------------------------------
2076   // Select active upper sectors
2077   //--------------------------------
2078
2079   //-----------------------------------------------------------------
2080   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2081   //-----------------------------------------------------------------
2082
2083   fSecUps[0] = s1;
2084   fSecUps[1] = s2;
2085   fSecUps[2] = s3;
2086   fSecUps[3] = s4;
2087   fSecUps[4] = s5;
2088   fSecUps[5] = s6;
2089   fSecUps[6] = s7;
2090   fSecUps[7] = s8;
2091   fSecUps[8] = s9;
2092   fSecUps[9] = s10;
2093   fSecUps[10] = s11;
2094   fSecUps[11] = s12;
2095 }
2096
2097 //_____________________________________________________________________________
2098 void AliTPC::SetSens(Int_t sens)
2099 {
2100
2101   //-------------------------------------------------------------
2102   // Activates/deactivates the sensitive strips at the center of
2103   // the pad row -- this is for the space-point resolution calculations
2104   //-------------------------------------------------------------
2105
2106   //-----------------------------------------------------------------
2107   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2108   //-----------------------------------------------------------------
2109
2110   fSens = sens;
2111 }
2112
2113 //_____________________________________________________________________________
2114 void AliTPC::Streamer(TBuffer &R__b)
2115 {
2116   //
2117   // Stream an object of class AliTPC.
2118   //
2119    if (R__b.IsReading()) {
2120       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
2121       AliDetector::Streamer(R__b);
2122       if (R__v < 2) return;
2123       R__b >> fNsectors;
2124       R__b >> fNclusters;
2125       R__b >> fNtracks;
2126       fClustersIndex = new Int_t[fNsectors+1];
2127       fDigitsIndex   = new Int_t[fNsectors+1];
2128    } else {
2129       R__b.WriteVersion(AliTPC::IsA());
2130       AliDetector::Streamer(R__b);
2131       R__b << fNsectors;
2132       R__b << fNclusters;
2133       R__b << fNtracks;
2134    }
2135 }
2136
2137 ClassImp(AliTPCcluster)
2138
2139 //_____________________________________________________________________________
2140 AliTPCcluster::AliTPCcluster(Float_t *hits, Int_t *lab)
2141 {
2142   //
2143   // Creates a simulated cluster for the TPC
2144   //
2145   fTracks[0]  = lab[0];
2146   fTracks[1]  = lab[1];
2147   fTracks[2]  = lab[2];
2148   fSector     = lab[3];
2149   fPadRow     = lab[4];
2150   fY          = hits[0];
2151   fZ          = hits[1];
2152   fQ          = hits[2];
2153   fSigmaY2    = hits[3];
2154   fSigmaZ2    = hits[4];
2155 }
2156
2157 //_____________________________________________________________________________
2158 void AliTPCcluster::GetXYZ(Float_t *x, const AliTPCParam *par) const
2159 {
2160   //
2161   // Transformation from local to global coordinate system
2162   //
2163   x[0]=par->GetPadRowRadii(fSector,fPadRow-1);
2164   x[1]=fY;
2165   x[2]=fZ;
2166   par->CRXYZtoXYZ(x,fSector,fPadRow-1,1);
2167 }
2168
2169 //_____________________________________________________________________________
2170 Int_t AliTPCcluster::Compare(TObject * o)
2171 {
2172   //
2173   // compare two clusters according y coordinata
2174   AliTPCcluster *cl= (AliTPCcluster *)o;
2175   if (fY<cl->fY) return -1;
2176   if (fY==cl->fY) return 0;
2177   return 1;
2178 }
2179
2180 Bool_t AliTPCcluster::IsSortable() const
2181 {
2182   //
2183   //make AliTPCcluster sortabale
2184   return kTRUE;
2185 }
2186
2187
2188
2189 ClassImp(AliTPCdigit)
2190
2191 //_____________________________________________________________________________
2192 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
2193   AliDigit(tracks)
2194 {
2195   //
2196   // Creates a TPC digit object
2197   //
2198   fSector     = digits[0];
2199   fPadRow     = digits[1];
2200   fPad        = digits[2];
2201   fTime       = digits[3];
2202   fSignal     = digits[4];
2203 }
2204
2205
2206 ClassImp(AliTPChit)
2207
2208 //_____________________________________________________________________________
2209 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
2210 AliHit(shunt,track)
2211 {
2212   //
2213   // Creates a TPC hit object
2214   //
2215   fSector     = vol[0];
2216   fPadRow     = vol[1];
2217   fX          = hits[0];
2218   fY          = hits[1];
2219   fZ          = hits[2];
2220   fQ          = hits[3];
2221 }
2222
2223
2224 ClassImp(AliTPCtrack)
2225
2226 //_____________________________________________________________________________
2227 AliTPCtrack::AliTPCtrack(Float_t *hits)
2228 {
2229   //
2230   // Default creator for a TPC reconstructed track object
2231   //
2232   ref=hits[0]; // This is dummy code !
2233 }
2234
2235 AliTPCtrack::AliTPCtrack(const AliTPCcluster& c,const TVector& xx,
2236                          const TMatrix& CC, const AliTPCParam *p):
2237   x(xx),C(CC),clusters(MAX_CLUSTER)
2238 {
2239   //
2240   // Standard creator for a TPC reconstructed track object
2241   //
2242   chi2=0.;
2243   int sec=c.fSector-1, row=c.fPadRow-1;
2244   ref=p->GetPadRowRadii(sec+1,row);
2245
2246   if (sec<24) {
2247     fAlpha=(sec%12)*alpha_low;
2248   } else {
2249     fAlpha=(sec%24)*alpha_up;
2250   }
2251   clusters.AddLast((AliTPCcluster*)(&c));
2252 }
2253
2254 //_____________________________________________________________________________
2255 AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : x(t.x), C(t.C),
2256   clusters(t.clusters.GetEntriesFast())
2257 {
2258   //
2259   // Copy creator for a TPC reconstructed track
2260   //
2261   ref=t.ref;
2262   chi2=t.chi2;
2263   fAlpha=t.fAlpha;
2264   int n=t.clusters.GetEntriesFast();
2265   for (int i=0; i<n; i++) clusters.AddLast(t.clusters.UncheckedAt(i));
2266 }
2267
2268 //_____________________________________________________________________________
2269 Double_t AliTPCtrack::GetY(Double_t xk) const
2270 {
2271   //
2272   //
2273   //
2274   Double_t c2=x(2)*xk - x(3);
2275   if (TMath::Abs(c2) >= 0.999) {
2276     if (*this>10) cerr<<*this<<" AliTPCtrack warning: No y for this x !\n";
2277     return 0.;
2278   }
2279   Double_t c1=x(2)*ref - x(3);
2280   Double_t r1=sqrt(1.-c1*c1), r2=sqrt(1.-c2*c2);
2281   Double_t dx=xk-ref;
2282   return x(0) + dx*(c1+c2)/(r1+r2);
2283 }
2284
2285 //_____________________________________________________________________________
2286 int AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
2287 {
2288   //
2289   // Propagate a TPC reconstructed track
2290   //
2291   if (TMath::Abs(x(2)*xk - x(3)) >= 0.999) {
2292     if (*this>3) cerr<<*this<<" AliTPCtrack warning: Propagation failed !\n";
2293     return 0;
2294   }
2295
2296   Double_t x1=ref, x2=x1+0.5*(xk-x1), dx=x2-x1, y1=x(0), z1=x(1);
2297   Double_t c1=x(2)*x1 - x(3), r1=sqrt(1.- c1*c1);
2298   Double_t c2=x(2)*x2 - x(3), r2=sqrt(1.- c2*c2);
2299
2300   x(0) += dx*(c1+c2)/(r1+r2);
2301   x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
2302
2303   TMatrix F(5,5); F.UnitMatrix();
2304   Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2;
2305   F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
2306   F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
2307   Double_t cr=c1*r2+c2*r1;
2308   F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
2309   F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
2310   F(1,4)= dx*cc/cr;
2311   TMatrix tmp(F,TMatrix::kMult,C);
2312   C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2313
2314   ref=x2;
2315
2316   //Multiple scattering******************
2317   Double_t ey=x(2)*ref - x(3);
2318   Double_t ex=sqrt(1-ey*ey);
2319   Double_t ez=x(4);
2320   TMatrix Q(5,5); Q=0.;
2321   Q(2,2)=ez*ez+ey*ey;   Q(2,3)=-ex*ey;       Q(2,4)=-ex*ez;
2322   Q(3,2)=Q(2,3);        Q(3,3)= ez*ez+ex*ex; Q(3,4)=-ey*ez;
2323   Q(4,2)=Q(2,4);        Q(4,3)= Q(3,4);      Q(4,4)=1.;
2324
2325   F=0;
2326   F(2,2)=-x(2)*ex;          F(2,3)=-x(2)*ey;
2327   F(3,2)=-ex*(x(2)*ref-ey); F(3,3)=-(1.+ x(2)*ref*ey - ey*ey);
2328   F(4,2)=-ez*ex;            F(4,3)=-ez*ey;           F(4,4)=1.;
2329
2330   tmp.Mult(F,Q);
2331   Q.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2332
2333   Double_t p2=GetPt()*GetPt()*(1.+x(4)*x(4));
2334   Double_t beta2=p2/(p2 + pm*pm);
2335   Double_t d=sqrt((x1-ref)*(x1-ref)+(y1-x(0))*(y1-x(0))+(z1-x(1))*(z1-x(1)));
2336   d*=2.;
2337   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
2338   Q*=theta2;
2339   C+=Q;
2340
2341   //Energy losses************************
2342   Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho;
2343   if (x1 < x2) dE=-dE;
2344   x(2)*=(1.- sqrt(p2+pm*pm)/p2*dE);
2345   //x(3)*=(1.- sqrt(p2+pm*pm)/p2*dE);
2346
2347   x1=ref; x2=xk; y1=x(0); z1=x(1);
2348   c1=x(2)*x1 - x(3); r1=sqrt(1.- c1*c1);
2349   c2=x(2)*x2 - x(3); r2=sqrt(1.- c2*c2);
2350
2351   x(0) += dx*(c1+c2)/(r1+r2);
2352   x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
2353
2354   F.UnitMatrix();
2355   rr=r1+r2; cc=c1+c2; xx=x1+x2;
2356   F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
2357   F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
2358   cr=c1*r2+c2*r1;
2359   F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
2360   F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
2361   F(1,4)= dx*cc/cr;
2362   tmp.Mult(F,C);
2363   C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2364
2365   ref=x2;
2366
2367   return 1;
2368 }
2369
2370 //_____________________________________________________________________________
2371 void AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm)
2372 {
2373   //
2374   // Propagate a reconstructed track from the vertex
2375   //
2376   Double_t c=x(2)*ref - x(3);
2377   Double_t tgf=-x(3)/(x(2)*x(0) + sqrt(1-c*c));
2378   Double_t snf=tgf/sqrt(1.+ tgf*tgf);
2379   Double_t xv=(x(3)+snf)/x(2);
2380   PropagateTo(xv,x0,rho,pm);
2381 }
2382
2383 //_____________________________________________________________________________
2384 void AliTPCtrack::Update(const AliTPCcluster *c, Double_t chisq)
2385 {
2386   //
2387   // Update statistics for a reconstructed TPC track
2388   //
2389   TMatrix H(2,5); H.UnitMatrix();
2390   TMatrix Ht(TMatrix::kTransposed,H);
2391   TVector m(2);   m(0)=c->fY; m(1)=c->fZ;
2392   TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
2393
2394   TMatrix tmp(H,TMatrix::kMult,C);
2395   TMatrix R(tmp,TMatrix::kMult,Ht); R+=V;
2396
2397   Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
2398   R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
2399   R(1,0)*=-1; R(0,1)=R(1,0);
2400   R*=1./det;
2401
2402   //R.Invert();
2403
2404   TMatrix K(C,TMatrix::kMult,Ht); K*=R;
2405
2406   TVector savex=x;
2407   x*=H; x-=m; x*=-1; x*=K; x+=savex;
2408   if (TMath::Abs(x(2)*ref-x(3)) >= 0.999) {
2409     if (*this>3) cerr<<*this<<" AliTPCtrack warning: Filtering failed !\n";
2410     x=savex;
2411     return;
2412   }
2413
2414   TMatrix saveC=C;
2415   C.Mult(K,tmp); C-=saveC; C*=-1;
2416
2417   clusters.AddLast((AliTPCcluster*)c);
2418   chi2 += chisq;
2419 }
2420
2421 //_____________________________________________________________________________
2422 int AliTPCtrack::Rotate(Double_t alpha)
2423 {
2424   //
2425   // Rotate a reconstructed TPC track
2426   //
2427   fAlpha += alpha;
2428
2429   Double_t x1=ref, y1=x(0);
2430   Double_t ca=cos(alpha), sa=sin(alpha);
2431   Double_t r1=x(2)*ref - x(3);
2432
2433   ref = x1*ca + y1*sa;
2434   x(0)=-x1*sa + y1*ca;
2435   x(3)=x(3)*ca + (x(2)*y1 + sqrt(1.- r1*r1))*sa;
2436
2437   Double_t r2=x(2)*ref - x(3);
2438   if (TMath::Abs(r2) >= 0.999) {
2439     if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !\n";
2440     return 0;
2441   }
2442
2443   Double_t y0=x(0) + sqrt(1.- r2*r2)/x(2);
2444   if ((x(0)-y0)*x(2) >= 0.) {
2445     if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !!!\n";
2446     return 0;
2447   }
2448
2449   TMatrix F(5,5); F.UnitMatrix();
2450   F(0,0)=ca;
2451   F(3,0)=x(2)*sa;
2452   F(3,2)=(y1 - r1*x1/sqrt(1.- r1*r1))*sa;
2453   F(3,3)= ca + sa*r1/sqrt(1.- r1*r1);
2454   TMatrix tmp(F,TMatrix::kMult,C);
2455   // Double_t dy2=C(0,0);
2456   C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
2457   // C(0,0)+=dy2*sa*sa*r1*r1/(1.- r1*r1);
2458   // C(1,1)+=dy2*sa*sa*x(4)*x(4)/(1.- r1*r1);
2459
2460   return 1;
2461 }
2462
2463 //_____________________________________________________________________________
2464 void AliTPCtrack::UseClusters() const
2465 {
2466   //
2467   //
2468   //
2469   int num_of_clusters=clusters.GetEntriesFast();
2470   for (int i=0; i<num_of_clusters; i++) {
2471     //if (i<=14) continue;
2472     AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i);
2473     c->Use();
2474   }
2475 }
2476
2477 //_____________________________________________________________________________
2478 Double_t AliTPCtrack::GetPredictedChi2(const AliTPCcluster *c) const
2479 {
2480   //
2481   // Calculate chi2 for a reconstructed TPC track
2482   //
2483   TMatrix H(2,5); H.UnitMatrix();
2484   TVector m(2);   m(0)=c->fY; m(1)=c->fZ;
2485   TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
2486   TVector res=x;  res*=H; res-=m; //res*=-1;
2487   TMatrix tmp(H,TMatrix::kMult,C);
2488   TMatrix R(tmp,TMatrix::kMult,TMatrix(TMatrix::kTransposed,H)); R+=V;
2489
2490   Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
2491   if (TMath::Abs(det) < 1.e-10) {
2492     if (*this>3) cerr<<*this<<" AliTPCtrack warning: Singular matrix !\n";
2493     return 1e10;
2494   }
2495   R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
2496   R(1,0)*=-1; R(0,1)=R(1,0);
2497   R*=1./det;
2498
2499   //R.Invert();
2500
2501   TVector r=res;
2502   res*=R;
2503   return r*res;
2504 }
2505
2506 //_____________________________________________________________________________
2507 int AliTPCtrack::GetLab() const
2508 {
2509   //
2510   //
2511   //
2512   int lab;
2513   struct {
2514     int lab;
2515     int max;
2516   } s[MAX_CLUSTER]={{0,0}};
2517
2518   int i;
2519   int num_of_clusters=clusters.GetEntriesFast();
2520   for (i=0; i<num_of_clusters; i++) {
2521     AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i);
2522     lab=TMath::Abs(c->fTracks[0]);
2523     int j;
2524     for (j=0; j<MAX_CLUSTER; j++)
2525       if (s[j].lab==lab || s[j].max==0) break;
2526     s[j].lab=lab;
2527     s[j].max++;
2528   }
2529
2530   int max=0;
2531   for (i=0; i<num_of_clusters; i++)
2532     if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
2533
2534   for (i=0; i<num_of_clusters; i++) {
2535     AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i);
2536     if (TMath::Abs(c->fTracks[1]) == lab ||
2537         TMath::Abs(c->fTracks[2]) == lab ) max++;
2538   }
2539
2540   if (1.-float(max)/num_of_clusters > 0.10) return -lab;
2541
2542   if (num_of_clusters < 6) return lab;
2543
2544   max=0;
2545   for (i=1; i<=6; i++) {
2546     AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(num_of_clusters-i);
2547     if (lab == TMath::Abs(c->fTracks[0]) ||
2548         lab == TMath::Abs(c->fTracks[1]) ||
2549         lab == TMath::Abs(c->fTracks[2])) max++;
2550   }
2551   if (max<3) return -lab;
2552
2553   return lab;
2554 }
2555
2556 //_____________________________________________________________________________
2557 void AliTPCtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const
2558 {
2559   //
2560   // Get reconstructed TPC track momentum
2561   //
2562   Double_t pt=0.3*FIELD/TMath::Abs(x(2))/100; // GeV/c
2563   Double_t r=x(2)*ref-x(3);
2564   Double_t y0=x(0) + sqrt(1.- r*r)/x(2);
2565   px=-pt*(x(0)-y0)*x(2);    //cos(phi);
2566   py=-pt*(x(3)-ref*x(2));   //sin(phi);
2567   pz=pt*x(4);
2568   Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha);
2569   py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha);
2570   px=tmp;
2571 }
2572
2573 //_____________________________________________________________________________
2574 //
2575 //     Classes for internal tracking use
2576 //
2577
2578 //_____________________________________________________________________________
2579 void AliTPCRow::InsertCluster(const AliTPCcluster* c)
2580 {
2581   //
2582   // Insert a cluster in the list
2583   //
2584   if (num_of_clusters==MAX_CLUSTER_PER_ROW) {
2585     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
2586   }
2587   if (num_of_clusters==0) {clusters[num_of_clusters++]=c; return;}
2588   int i=Find(c->fY);
2589   memmove(clusters+i+1 ,clusters+i,(num_of_clusters-i)*sizeof(AliTPCcluster*));
2590   clusters[i]=c; num_of_clusters++;
2591 }
2592
2593 //_____________________________________________________________________________
2594 int AliTPCRow::Find(Double_t y) const
2595 {
2596   //
2597   //
2598   //
2599   if (y <= clusters[0]->fY) return 0;
2600   if (y > clusters[num_of_clusters-1]->fY) return num_of_clusters;
2601   int b=0, e=num_of_clusters-1, m=(b+e)/2;
2602   for (; b<e; m=(b+e)/2) {
2603     if (y > clusters[m]->fY) b=m+1;
2604     else e=m;
2605   }
2606   return m;
2607 }
2608