TPC/AliTPC.cxx

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