TPC/AliTPC.cxx

   1 ///////////////////////////////////////////////////////////////////////////////
   2 //                                                                           //
   3 //  Time Projection Chamber                                                  //
   4 //  This class contains the basic functions for the Time Projection Chamber  //
   5 //  detector. Functions specific to one particular geometry are              //
   6 //  contained in the derived classes                                         //
   7 //                                                                           //
   8 //Begin_Html
   9 /*
  10 <img src="gif/AliTPCClass.gif">
  11 */
  12 //End_Html
  13 //                                                                           //
  14 //                                                                           //
  15 ///////////////////////////////////////////////////////////////////////////////
  16
  17 #include <TMath.h>
  18 #include <TRandom.h>
  19 #include <TVector.h>
  20 #include <TGeometry.h>
  21 #include <TNode.h>
  22 #include <TTUBS.h>
  23 #include <TObjectTable.h>
  24 #include "GParticle.h"
  25 #include "AliTPC.h"
  26 #include "AliRun.h"
  27 #include <iostream.h>
  28 #include <fstream.h>
  29 #include "AliMC.h"
  30
  31 ClassImp(AliTPC)
  32
  33 //_____________________________________________________________________________
  34 AliTPC::AliTPC()
  35 {
  36   //
  37   // Default constructor
  38   //
  39   fIshunt   = 0;
  40   fClusters = 0;
  41   fHits     = 0;
  42   fDigits   = 0;
  43   fTracks   = 0;
  44   fNsectors = 0;
  45   fNtracks  = 0;
  46   fNclusters= 0;
  47 }
  48
  49 //_____________________________________________________________________________
  50 AliTPC::AliTPC(const char *name, const char *title)
  51       : AliDetector(name,title)
  52 {
  53   //
  54   // Standard constructor
  55   //
  56
  57   //
  58   // Initialise arrays of hits and digits
  59   fHits     = new TClonesArray("AliTPChit",  176);
  60   fDigits   = new TClonesArray("AliTPCdigit",10000);
  61   //
  62   // Initialise counters
  63   fClusters = 0;
  64   fTracks   = 0;
  65   fNsectors = 72;
  66   fNtracks  = 0;
  67   fNclusters= 0;
  68   fDigitsIndex = new Int_t[fNsectors+1];
  69   fClustersIndex = new Int_t[fNsectors+1];
  70   //
  71   fIshunt     =  0;
  72   //
  73   // Initialise color attributes
  74   SetMarkerColor(kYellow);
  75 }
  76
  77 //_____________________________________________________________________________
  78 AliTPC::~AliTPC()
  79 {
  80   //
  81   // TPC destructor
  82   //
  83   fIshunt   = 0;
  84   delete fHits;
  85   delete fDigits;
  86   delete fClusters;
  87   delete fTracks;
  88   if (fDigitsIndex)   delete [] fDigitsIndex;
  89   if (fClustersIndex) delete [] fClustersIndex;
  90 }
  91
  92 //_____________________________________________________________________________
  93 void AliTPC::AddCluster(Float_t *hits, Int_t *tracks)
  94 {
  95   //
  96   // Add a simulated cluster to the list
  97   //
  98   if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000);
  99   TClonesArray &lclusters = *fClusters;
 100   new(lclusters[fNclusters++]) AliTPCcluster(hits,tracks);
 101 }
 102
 103 //_____________________________________________________________________________
 104 void AliTPC::AddCluster(const AliTPCcluster &c)
 105 {
 106   //
 107   // Add a simulated cluster copy to the list
 108   //
 109   if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000);
 110   TClonesArray &lclusters = *fClusters;
 111   new(lclusters[fNclusters++]) AliTPCcluster(c);
 112 }
 113
 114 //_____________________________________________________________________________
 115 void AliTPC::AddDigit(Int_t *tracks, Int_t *digits)
 116 {
 117   //
 118   // Add a TPC digit to the list
 119   //
 120   TClonesArray &ldigits = *fDigits;
 121   new(ldigits[fNdigits++]) AliTPCdigit(tracks,digits);
 122 }
 123
 124 //_____________________________________________________________________________
 125 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 126 {
 127   //
 128   // Add a hit to the list
 129   //
 130   TClonesArray &lhits = *fHits;
 131   new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
 132 }
 133
 134 //_____________________________________________________________________________
 135 void AliTPC::AddTrack(Float_t *hits)
 136 {
 137   //
 138   // Add a track to the list of tracks
 139   //
 140   TClonesArray &ltracks = *fTracks;
 141   new(ltracks[fNtracks++]) AliTPCtrack(hits);
 142 }
 143
 144 //_____________________________________________________________________________
 145 void AliTPC::AddTrack(const AliTPCtrack& t)
 146 {
 147   //
 148   // Add a track copy to the list of tracks
 149   //
 150   if(!fTracks) fTracks=new TClonesArray("AliTPCtrack",10000);
 151   TClonesArray &ltracks = *fTracks;
 152   new(ltracks[fNtracks++]) AliTPCtrack(t);
 153 }
 154
 155 //_____________________________________________________________________________
 156 void AliTPC::BuildGeometry()
 157 {
 158   //
 159   // Build TPC ROOT TNode geometry for the event display
 160   //
 161   TNode *Node, *Top;
 162   TTUBS *tubs;
 163   Int_t i;
 164   const int kColorTPC=19;
 165   char name[5], title[18];
 166   const Double_t kDegrad=TMath::Pi()/180;
 167   const Double_t loAng=30;
 168   const Double_t hiAng=15;
 169   const Int_t nLo = Int_t (360/loAng+0.5);
 170   const Int_t nHi = Int_t (360/hiAng+0.5);
 171   const Double_t loCorr = 1/TMath::Cos(0.5*loAng*kDegrad);
 172   const Double_t hiCorr = 1/TMath::Cos(0.5*hiAng*kDegrad);
 173   //
 174   // Get ALICE top node
 175   Top=gAlice->GetGeometry()->GetNode("alice");
 176   //
 177   // Inner sectors
 178   for(i=0;i<nLo;i++) {
 179     sprintf(name,"LS%2.2d",i);
 180     sprintf(title,"TPC low sector %d",i);
 181     tubs = new TTUBS(name,title,"void",88*loCorr,136*loCorr,250,loAng*(i-0.5),loAng*(i+0.5));
 182     tubs->SetNumberOfDivisions(1);
 183     Top->cd();
 184     Node = new TNode(name,title,name,0,0,0,"");
 185     Node->SetLineColor(kColorTPC);
 186     fNodes->Add(Node);
 187   }
 188   // Outer sectors
 189   for(i=0;i<nHi;i++) {
 190     sprintf(name,"US%2.2d",i);
 191     sprintf(title,"TPC upper sector %d",i);
 192     tubs = new TTUBS(name,title,"void",142*hiCorr,250*hiCorr,250,hiAng*(i-0.5),hiAng*(i+0.5));
 193     tubs->SetNumberOfDivisions(1);
 194     Top->cd();
 195     Node = new TNode(name,title,name,0,0,0,"");
 196     Node->SetLineColor(kColorTPC);
 197     fNodes->Add(Node);
 198   }
 199 }
 200
 201
 202 //_____________________________________________________________________________
 203 void AliTPC::CreateList(Int_t *tracks,Float_t signal[][MAXTPCTBK+1],
 204                         Int_t ntr,Int_t time)
 205 {
 206   //
 207   // Creates list of tracks contributing to a given digit
 208   // Only the 3 most significant tracks are taken into account
 209   //
 210
 211   Int_t i,j;
 212
 213   for(i=0;i<3;i++) tracks[i]=-1;
 214
 215   //
 216   //  Loop over signals, only 3 times
 217   //
 218
 219   Float_t qmax;
 220   Int_t jmax;
 221   Int_t jout[3] = {-1,-1,-1};
 222
 223   for(i=0;i<3;i++){
 224     qmax=0.;
 225     jmax=0;
 226
 227     for(j=0;j<ntr;j++){
 228
 229       if((i == 1 && j == jout[i-1])
 230          ||(i == 2 && (j == jout[i-1] || j == jout[i-2]))) continue;
 231
 232       if(signal[j][time] > qmax) {
 233         qmax = signal[j][time];
 234         jmax=j;
 235       }
 236     }
 237
 238     if(qmax > 0.) {
 239       tracks[i]=jmax;
 240       jout[i]=jmax;
 241     }
 242
 243   }
 244
 245   for(i=0;i<3;i++){
 246     if(tracks[i] < 0){
 247       tracks[i]=0;
 248     }
 249     else {
 250       tracks[i]=(Int_t) signal[tracks[i]][0]; // first element is a track number
 251     }
 252   }
 253 }
 254
 255 //_____________________________________________________________________________
 256 void AliTPC::DigSignal(Int_t isec,Int_t irow,TObjArray *pointer)
 257 {
 258   //
 259   // Digitalise TPC signal
 260   //
 261   Int_t pad_c,n_of_pads;
 262   Int_t pad_number;
 263
 264   n_of_pads = (isec < 25) ? npads_low[irow] : npads_up[irow];
 265   pad_c=(n_of_pads+1)/2; // this is the "central" pad for a row
 266
 267   Int_t track,idx;
 268   Int_t entries;
 269   TVector *pp;
 270   TVector *ppp;
 271   Float_t y,yy,z;
 272   Float_t pad_signal = 0;
 273   Float_t signal[MAXTPCTBK]; // Integrated signal over all tracks
 274   Float_t signal_tr[100][MAXTPCTBK+1]; // contribution from less than 50 tracks
 275   Int_t flag; // flag indicating a track contributing to a pad signal
 276
 277   //
 278
 279   Int_t ntracks = pointer->GetEntriesFast();
 280
 281   if(ntracks == 0) return; // no signal on this row!
 282
 283   //--------------------------------------------------------------
 284   // For each electron calculate the pad number and the avalanche
 285   //             This is only once per pad row
 286   //--------------------------------------------------------------
 287
 288   TVector **el = new TVector* [ntracks]; // each track is a new vector
 289
 290   TObjArray *arr = new TObjArray; // array of tracks for this row
 291
 292   for(track=0;track<ntracks;track++) {
 293     pp = (TVector*) pointer->At(track);
 294     entries = pp->GetNrows();
 295     el[track] = new TVector(entries-1);
 296     TVector &v1 = *el[track];
 297     TVector &v2 = *pp;
 298
 299     for(idx=0;idx<entries-2;idx+=2)
 300       {
 301         y=v2(idx+1);
 302         yy=TMath::Abs(y);
 303
 304         Float_t range=((n_of_pads-1)/2 + 0.5)*pad_pitch_w;
 305         //
 306         // Pad number and pad range
 307         //
 308         if(yy < 0.5*pad_pitch_w){
 309           pad_number=pad_c;
 310         }
 311         else if (yy < range){
 312           pad_number=(Int_t) ((yy-0.5*pad_pitch_w)/pad_pitch_w +1.);
 313           pad_number=(Int_t) (pad_number*TMath::Sign(1.,(double) y)+pad_c);
 314         }
 315         else{
 316           pad_number=0;
 317         }
 318
 319         v1(idx) = (Float_t) pad_number;
 320
 321         // Avalanche, taking the fluctuations into account
 322
 323         Int_t gain_fluct = (Int_t) (-gas_gain*TMath::Log(gRandom->Rndm()));
 324         v1(idx+1)= (Float_t) gain_fluct;
 325
 326       } // end of loop over electrons
 327
 328     arr->Add(el[track]); // add the vector to the array
 329
 330   }// end of loop over tracks
 331
 332   delete [] el; //  delete an array of pointers
 333
 334   //-------------------------------------------------------------
 335   //  Calculation of signal for every pad
 336   //-------------------------------------------------------------
 337
 338   //-------------------------------------------------------------
 339   // Loop over pads
 340   //-------------------------------------------------------------
 341
 342
 343   for(Int_t np=1;np<n_of_pads+1;np++)
 344     {
 345       for(Int_t l =0;l<MAXTPCTBK;l++) signal[l]=0.; // set signals for this pad to 0
 346
 347       for(Int_t k1=0;k1<100;k1++){
 348         for(Int_t k2=0;k2<MAXTPCTBK+1;k2++) signal_tr[k1][k2]=0.;
 349       }
 350       Int_t track_counter=0;
 351       //
 352
 353       //---------------------------------------------------------
 354       // Loop over tracks
 355       // --------------------------------------------------------
 356
 357       for(track=0;track<ntracks;track++)
 358         {
 359           flag = 0;
 360           pp = (TVector*) pointer->At(track);
 361           ppp = (TVector*) arr->At(track);
 362
 363           TVector &v1 = *pp;
 364           TVector &v2 = *ppp;
 365
 366           entries = pp->GetNrows();
 367
 368
 369           //----------------------------------------------------
 370           // Loop over electrons
 371           //----------------------------------------------------
 372
 373           for(idx=0;idx<entries-2;idx+=2)
 374             {
 375
 376               pad_number = (Int_t) v2(idx);
 377
 378               if(pad_number == 0) continue; // skip electrons outside range
 379
 380               Int_t pad_dist = pad_number-np;
 381               Int_t abs_dist = TMath::Abs(pad_dist);
 382
 383               if(abs_dist > 3) continue; // beyond signal range
 384
 385               y=  v1(idx+1);
 386               z = v1(idx+2);
 387
 388               Float_t dist = y-(pad_number-pad_c)*pad_pitch_w;
 389
 390               //----------------------------------------------
 391               // Calculate the signal induced on a pad "np"
 392               //----------------------------------------------
 393
 394               if(pad_dist < 0) dist = -dist;
 395
 396               switch((Int_t) abs_dist){
 397               case 0 : pad_signal = P4(dist); // electron is on pad "np"
 398                           break;
 399               case 1 : pad_signal = P3(dist); // electron is 1 pad away
 400                 break;
 401               case 2 : pad_signal = P2(dist); // electron is 2 pads away
 402                 break;
 403               case 3 : pad_signal = P1(dist); // electron is 3 pads away
 404               }
 405
 406               //---------------------------------
 407               //  Multiply by a gas gain
 408               //---------------------------------
 409
 410               pad_signal=pad_signal*v2(idx+1);
 411
 412               flag = 1;
 413
 414
 415               //-----------------------------------------------
 416               //  Sample the pad signal in time
 417               //-----------------------------------------------
 418
 419               Float_t t_drift = (z_end-TMath::Abs(z))/v_drift; // drift time
 420
 421               Float_t t_offset = t_drift-t_sample*(Int_t)(t_drift/t_sample);
 422               Int_t first_bucket = (Int_t) (t_drift/t_sample+1.);
 423
 424               for(Int_t bucket = 1;bucket<6;bucket++){
 425                 Float_t time = (bucket-1)*t_sample+t_offset;
 426                 Int_t time_idx = first_bucket+bucket-1;
 427                 Float_t ampl = pad_signal*TimeRes(time);
 428                 if (time_idx > MAXTPCTBK) break; //Y.Belikov
 429                 if (track_counter >=100) break;  //Y.Belikov
 430
 431                 signal_tr[track_counter][time_idx] += ampl; // single track only
 432                 signal[time_idx-1] += ampl;  // fill a signal array for this pad
 433
 434               } // end of time sampling
 435
 436             } // end of loop over electrons for a current track
 437
 438           //-----------------------------------------------
 439           //  add the track number and update the counter
 440           //  if it gave a nonzero contribution to the pad
 441           //-----------------------------------------------
 442           if(flag != 0 && track_counter < 100){
 443           signal_tr[track_counter][0] = v1(0);
 444           track_counter++;      // counter is looking at the NEXT track!
 445           }
 446
 447         } // end of loop over tracks
 448
 449       //----------------------------------------------
 450       //  Fill the Digits for this pad
 451       //----------------------------------------------
 452
 453       Int_t tracks[3];
 454       Int_t digits[5];
 455
 456       digits[0]=isec; // sector number
 457       digits[1]=irow+1; // row number
 458       digits[2]=np;   // pad number
 459
 460       Float_t q;
 461
 462       for(Int_t time = 0;time<MAXTPCTBK;time++){
 463         digits[3] = time+1; // time bucket
 464
 465         q = signal[time];
 466         q = gRandom->Gaus(q,sigma_noise); // apply noise
 467
 468         q *= (q_el*1.e15); // convert to fC
 469         q *= chip_gain; // convert to mV
 470         q *= (adc_sat/dyn_range); // convert to ADC counts
 471
 472         if(q < zero_sup) continue; // do not fill "zeros"
 473         if(q > adc_sat) q = adc_sat;  // saturation
 474         digits[4] = (Int_t) q;                   // ADC counts
 475
 476         //--------------------------------------------------
 477         //    "Real signal" or electronics noise
 478         //--------------------------------------------------
 479
 480         if(signal[time] > 0.){
 481
 482           //-----------------------------------------------------
 483           //  Create a list of tracks contributing to this digit
 484           //  If the digit results from a noise, track number is 0
 485           //-----------------------------------------------------
 486
 487           CreateList(tracks,signal_tr,track_counter,time);
 488         }
 489         else {
 490           for(Int_t ii=0;ii<3;ii++) tracks[ii]=0;
 491         }
 492
 493         AddDigit(tracks,digits);
 494
 495
 496       } // end of digits for this pad
 497
 498     } // end of loop over pads
 499
 500   arr->Delete(); // delete objects in this array
 501
 502   delete arr; // deletes the TObjArray itselves
 503
 504 }
 505
 506 //_____________________________________________________________________________
 507 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
 508 {
 509   //
 510   // Calculate distance from TPC to mouse on the display
 511   // Dummy procedure
 512   //
 513   return 9999;
 514 }
 515
 516 //_____________________________________________________________________________
 517 const int MAX_CLUSTER=nrow_low+nrow_up;
 518 const int S_MAXSEC=24;
 519 const int L_MAXSEC=48;
 520 const int ROWS_TO_SKIP=21;
 521 const Float_t MAX_CHI2=15.;
 522 const Float_t THRESHOLD=8*zero_sup;
 523
 524 //_____________________________________________________________________________
 525 static Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
 526 {
 527   //
 528   // Calculate spread in Y
 529   //
 530   pt=TMath::Abs(pt)*1000.;
 531   Double_t x=r/pt;
 532   tgl=TMath::Abs(tgl);
 533   Double_t s=a_rphi - b_rphi*r*tgl + c_rphi*x*x + d_rphi*x;
 534   if (s<0.4e-3) s=0.4e-3;
 535   return s;
 536 }
 537
 538 //_____________________________________________________________________________
 539 static Double_t SigmaZ2(Double_t r, Double_t tgl)
 540 {
 541   //
 542   // Calculate spread in Z
 543   //
 544   tgl=TMath::Abs(tgl);
 545   Double_t s=a_z - b_z*r*tgl + c_z*tgl*tgl;
 546   if (s<0.4e-3) s=0.4e-3;
 547   return s;
 548 }
 549
 550 //_____________________________________________________________________________
 551 inline Double_t f1(Double_t x1,Double_t y1,   //C
 552                    Double_t x2,Double_t y2,
 553                    Double_t x3,Double_t y3)
 554 {
 555   //
 556   // Function f1
 557   //
 558   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 559   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 560                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 561   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 562                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 563   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 564
 565   return -xr*yr/sqrt(xr*xr+yr*yr);
 566 }
 567
 568
 569 //_____________________________________________________________________________
 570 inline Double_t f2(Double_t x1,Double_t y1,  //eta=C*x0
 571                    Double_t x2,Double_t y2,
 572                    Double_t x3,Double_t y3)
 573 {
 574   //
 575   // Function f2
 576   //
 577   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 578   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 579                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 580   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 581                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 582   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 583
 584   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 585 }
 586
 587 //_____________________________________________________________________________
 588 inline Double_t f3(Double_t x1,Double_t y1,  //tgl
 589                    Double_t x2,Double_t y2,
 590                    Double_t z1,Double_t z2)
 591 {
 592   //
 593   // Function f3
 594   //
 595   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 596 }
 597
 598 //_____________________________________________________________________________
 599 static int FindProlongation(AliTPCtrack& t, const AliTPCSector *sec,
 600                             int s, int ri, int rf=0)
 601 {
 602   //
 603   // Propagate track
 604   //
 605   int try_again=ROWS_TO_SKIP;
 606   Double_t alpha=sec->GetAlpha();
 607   int ns=int(2*TMath::Pi()/alpha)+1;
 608   for (int nr=ri; nr>=rf; nr--) {
 609     Double_t x=sec[s].GetX(nr), ymax=sec[s].GetMaxY(nr);
 610     if (!t.PropagateTo(x)) break;
 611
 612     const AliTPCcluster *cl=0;
 613     Double_t max_chi2=MAX_CHI2;
 614     const AliTPCRow& row=sec[s][nr];
 615     Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
 616     Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
 617     Double_t road=3.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
 618
 619     if (road>30) {
 620       if (t>3) cerr<<t<<" AliTPCtrack warning: Too broad road !\n";
 621       break;
 622     }
 623
 624     if (row)
 625       for (int i=row.Find(y-road); i<row; i++) {
 626         AliTPCcluster* c=(AliTPCcluster*)(row[i]);
 627         if (c->fY > y+road) break;
 628         if (c->IsUsed()) continue;
 629         if ((c->fZ - z)*(c->fZ - z) > 9.*(t.GetSigmaZ2() + sz2)) continue;
 630         Double_t chi2=t.GetPredictedChi2(c);
 631         if (chi2 > max_chi2) continue;
 632         max_chi2=chi2;
 633         cl=c;
 634       }
 635     if (cl) {
 636       t.Update(cl,max_chi2);
 637       try_again=ROWS_TO_SKIP;
 638     } else {
 639       if (try_again--) {
 640         if (y > ymax) {
 641           s = (s+1) % ns;
 642           if (!t.Rotate(alpha)) break;
 643         } else
 644           if (y <-ymax) {
 645             s = (s-1+ns) % ns;
 646             if (!t.Rotate(-alpha)) break;
 647           };
 648         continue;
 649       }
 650       break;
 651     }
 652   }
 653   return s;
 654 }
 655
 656 //_____________________________________________________________________________
 657 static void MakeSeeds(TObjArray& seeds,const AliTPCSector* sec,int i1,int i2)
 658 {
 659   //
 660   // Find seed for tracking
 661   //
 662   TMatrix C(5,5); TVector x(5);
 663   int max_sec=L_MAXSEC/2;
 664   for (int ns=0; ns<max_sec; ns++) {
 665     int nl=sec[(ns-1+max_sec)%max_sec][i2];
 666     int nm=sec[ns][i2];
 667     int nu=sec[(ns+1)%max_sec][i2];
 668     Double_t alpha=sec[ns].GetAlpha();
 669     const AliTPCRow& r1=sec[ns][i1];
 670     for (int is=0; is < r1; is++) {
 671       Double_t x1=sec[ns].GetX(i1), y1=r1[is]->fY, z1=r1[is]->fZ;
 672       for (int js=0; js < nl+nm+nu; js++) {
 673         const AliTPCcluster *cl;
 674         Double_t cs,sn;
 675         //int ks;
 676
 677         if (js<nl) {
 678           //ks=(ns-1+max_sec)%max_sec;
 679           const AliTPCRow& r2=sec[(ns-1+max_sec)%max_sec][i2];
 680           cl=r2[js];
 681           cs=cos(alpha); sn=sin(alpha);
 682         } else
 683           if (js<nl+nm) {
 684             //ks=ns;
 685             const AliTPCRow& r2=sec[ns][i2];
 686             cl=r2[js-nl];
 687             cs=1; sn=0.;
 688           } else {
 689             //ks=(ns+1)%max_sec;
 690             const AliTPCRow& r2=sec[(ns+1)%max_sec][i2];
 691             cl=r2[js-nl-nm];
 692             cs=cos(alpha); sn=-sin(alpha);
 693           }
 694         Double_t x2=sec[ns].GetX(i2), y2=cl->fY, z2=cl->fZ;
 695         //if (z1*z2 < 0) continue;
 696         //if (TMath::Abs(z1) < TMath::Abs(z2)) continue;
 697
 698         Double_t tmp= x2*cs+y2*sn;
 699         y2 =-x2*sn+y2*cs;
 700         x2=tmp;
 701
 702         x(0)=y1;
 703         x(1)=z1;
 704         x(2)=f1(x1,y1,x2,y2,0.,0.);
 705         x(3)=f2(x1,y1,x2,y2,0.,0.);
 706         x(4)=f3(x1,y1,x2,y2,z1,z2);
 707
 708         if (TMath::Abs(x(2)*x1-x(3)) >= 0.999) continue;
 709
 710         if (TMath::Abs(x(4)) > 1.2) continue;
 711
 712         Double_t a=asin(x(3));
 713         /*
 714           Double_t tgl1=z1*x(2)/(a+asin(x(2)*x1-x(3)));
 715           Double_t tgl2=z2*x(2)/(a+asin(x(2)*x2-x(3)));
 716           Double_t ratio=2*(tgl1-tgl2)/(tgl1+tgl2);
 717           if (TMath::Abs(ratio)>0.0170) continue; //or > 0.005
 718         */
 719         Double_t zv=z1 - x(4)/x(2)*(a+asin(x(2)*x1-x(3)));
 720         if (TMath::Abs(zv)>33.) continue;
 721
 722
 723         TMatrix X(6,6); X=0.;
 724         X(0,0)=r1[is]->fSigmaY2; X(1,1)=r1[is]->fSigmaZ2;
 725         X(2,2)=cl->fSigmaY2;     X(3,3)=cl->fSigmaZ2;
 726         X(4,4)=3./12.; X(5,5)=3./12.;
 727         TMatrix F(5,6); F.UnitMatrix();
 728         Double_t sy=sqrt(X(0,0)), sz=sqrt(X(1,1));
 729         F(2,0)=(f1(x1,y1+sy,x2,y2,0.,0.)-x(2))/sy;
 730         F(2,2)=(f1(x1,y1,x2,y2+sy,0.,0.)-x(2))/sy;
 731         F(2,4)=(f1(x1,y1,x2,y2,0.,0.+sy)-x(2))/sy;
 732         F(3,0)=(f2(x1,y1+sy,x2,y2,0.,0.)-x(3))/sy;
 733         F(3,2)=(f2(x1,y1,x2,y2+sy,0.,0.)-x(3))/sy;
 734         F(3,4)=(f2(x1,y1,x2,y2,0.,0.+sy)-x(3))/sy;
 735         F(4,0)=(f3(x1,y1+sy,x2,y2,z1,z2)-x(4))/sy;
 736         F(4,1)=(f3(x1,y1,x2,y2,z1+sz,z2)-x(4))/sz;
 737         F(4,2)=(f3(x1,y1,x2,y2+sy,z1,z2)-x(4))/sy;
 738         F(4,3)=(f3(x1,y1,x2,y2,z1,z2+sz)-x(4))/sz;
 739         F(4,4)=0;
 740         F(3,3)=0;
 741
 742         TMatrix t(F,TMatrix::kMult,X);
 743         C.Mult(t,TMatrix(TMatrix::kTransposed,F));
 744
 745         TrackSeed *track=new TrackSeed(*(r1[is]),x,C);
 746         FindProlongation(*track,sec,ns,i1-1,i2);
 747         int ii=(i1-i2)/2;
 748         if (*track >= ii) {seeds.AddLast(track); continue;}
 749         else delete track;
 750       }
 751     }
 752   }
 753 }
 754
 755 //_____________________________________________________________________________
 756 void AliTPC::Clusters2Tracks()
 757 {
 758   //
 759   // TPC Track finder from clusters.
 760   //
 761   if (!fClusters) return;
 762   AliTPCSSector ssec[S_MAXSEC/2];
 763   AliTPCLSector lsec[L_MAXSEC/2];
 764   int ncl=fClusters->GetEntriesFast();
 765   while (ncl--) {
 766     AliTPCcluster *c=(AliTPCcluster*)fClusters->UncheckedAt(ncl);
 767     int sec=int(c->fSector), row=int(c->fPadRow);
 768
 769     if (sec<24) {
 770       if (row<0 || row>nrow_low) {cerr<<"low !!!"<<row<<endl; continue;}
 771       ssec[sec%12][row].InsertCluster(c);
 772     } else {
 773       if (row<0 || row>nrow_up ) {cerr<<"up  !!!"<<row<<endl; continue;}
 774       sec -= 24;
 775       lsec[sec%24][row].InsertCluster(c);
 776     }
 777   }
 778
 779
 780   TObjArray seeds(20000);
 781   MakeSeeds(seeds,lsec,nrow_up-1,nrow_up-1-8);
 782   MakeSeeds(seeds,lsec,nrow_up-1-4,nrow_up-1-4-8);
 783
 784   seeds.Sort();
 785
 786   int found=0;
 787   int nseed=seeds.GetEntriesFast();
 788
 789   for (int s=0; s<nseed; s++) {
 790     AliTPCtrack& t=*((AliTPCtrack*)seeds.UncheckedAt(s));
 791     Double_t alpha=t.GetAlpha();
 792     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 793     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 794     int ns=int(alpha/lsec->GetAlpha() + 0.5);
 795
 796     Double_t x=t.GetX();
 797     int nr;
 798     if (x<pad_row_up[nrow_up-1-4-7]) nr=nrow_up-1-4-8;
 799     else if (x<pad_row_up[nrow_up-1-7]) nr=nrow_up-1-8;
 800     else {cerr<<x<<" =x !!!\n"; continue;}
 801
 802     int ls=FindProlongation(t,lsec,ns,nr-1);
 803
 804     //   if (t < 25) continue;
 805
 806     x=t.GetX(); alpha=lsec[ls].GetAlpha();          //
 807     Double_t phi=ls*alpha + atan(t.GetY()/x);       // Find S-sector
 808     int ss=int(0.5*(phi/alpha+1));                  //
 809     alpha *= 2*(ss-0.5*ls);                         // and rotation angle
 810     if (!t.Rotate(alpha)) continue;                 //
 811     ss %= (S_MAXSEC/2);                             //
 812
 813     ss=FindProlongation(t,ssec,ss,nrow_low-1);
 814     if (t < 30) continue;
 815
 816     AddTrack(t);
 817     t.UseClusters();
 818     cerr<<found++<<'\r';
 819   }
 820 }
 821
 822 //_____________________________________________________________________________
 823 void AliTPC::CreateMaterials()
 824 {
 825   //
 826   // Create Materials for for TPC
 827   // Origin M.Kowalski
 828   //
 829
 830   AliMC* pMC = AliMC::GetMC();
 831
 832   Int_t ISXFLD=gAlice->Field()->Integ();
 833   Float_t SXMGMX=gAlice->Field()->Max();
 834
 835   Float_t absl, radl, a, d, z;
 836   Float_t dg;
 837   Float_t x0ne;
 838   Float_t buf[1];
 839   Int_t nbuf;
 840
 841   // --- Methane (CH4) ---
 842   Float_t am[2] = { 12.,1. };
 843   Float_t zm[2] = { 6.,1. };
 844   Float_t wm[2] = { 1.,4. };
 845   Float_t dm    = 7.17e-4;
 846   // --- The Neon CO2 90/10 mixture ---
 847   Float_t ag[2] = { 20.18 };
 848   Float_t zg[2] = { 10. };
 849   Float_t wg[2] = { .8,.2 };
 850   Float_t dne   = 9e-4;        // --- Neon density in g/cm3 ---
 851   // --- Mylar (C5H4O2) ---
 852   Float_t amy[3] = { 12.,1.,16. };
 853   Float_t zmy[3] = { 6.,1.,8. };
 854   Float_t wmy[3] = { 5.,4.,2. };
 855   Float_t dmy    = 1.39;
 856   // --- CO2 ---
 857   Float_t ac[2] = { 12.,16. };
 858   Float_t zc[2] = { 6.,8. };
 859   Float_t wc[2] = { 1.,2. };
 860   Float_t dc    = .001977;
 861   // --- Carbon density and radiation length ---
 862   Float_t densc = 2.265;
 863   Float_t radlc = 18.8;
 864   // --- Silicon ---
 865   Float_t asi   = 28.09;
 866   Float_t zsi   = 14.;
 867   Float_t desi  = 2.33;
 868   Float_t radsi = 9.36;
 869
 870   // --- Define the various materials for GEANT ---
 871   AliMaterial(0, "Al $", 26.98, 13., 2.7, 8.9, 37.2);
 872   x0ne = 28.94 / dne;
 873   AliMaterial(1, "Ne $", 20.18, 10., dne, x0ne, 999.);
 874
 875   // --  Methane, defined by the proportions of atoms
 876
 877   AliMixture(2, "Methane$", am, zm, dm, -2, wm);
 878
 879   // --- CO2, defined by the proportion of atoms
 880
 881   AliMixture(7, "CO2$", ac, zc, dc, -2, wc);
 882
 883   // --  Get A,Z etc. for CO2
 884
 885   char namate[21];
 886   pMC->Gfmate((*fIdmate)[7], namate, a, z, d, radl, absl, buf, nbuf);
 887
 888   ag[1] = a;
 889   zg[1] = z;
 890   dg = dne * .9 + dc * .1;
 891
 892   //-------------------------------------
 893
 894   // --  Create Ne/CO2 90/10 mixture
 895
 896   AliMixture(3, "Gas-mixt $", ag, zg, dg, 2, wg);
 897   AliMixture(4, "Gas-mixt $", ag, zg, dg, 2, wg);
 898
 899   AliMaterial(5, "G10$", 20., 10., 1.7, 19.4, 999.);
 900   AliMixture(6, "Mylar$", amy, zmy, dmy, -3, wmy);
 901
 902   a = ac[0];
 903   z = zc[0];
 904   AliMaterial(8, "Carbon", a, z, densc, radlc, 999.);
 905
 906   AliMaterial(9, "Silicon", asi, zsi, desi, radsi, 999.);
 907   AliMaterial(99, "Air$", 14.61, 7.3, .001205, 30420., 67500.);
 908
 909   AliMedium(400, "Al wall$",  0, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1);
 910   AliMedium(402, "Gas mix1$", 3, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
 911   AliMedium(403, "Gas mix2$", 3, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
 912   AliMedium(404, "Gas mix3$", 4, 1, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
 913   AliMedium(405, "G10 pln$",  5, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1 );
 914   AliMedium(406, "Mylar  $",  6, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
 915   AliMedium(407, "CO2    $",  7, 0, ISXFLD, SXMGMX, 10., .01,.1, .01,  .01);
 916   AliMedium(408, "Carbon $",  8, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1 );
 917   AliMedium(409, "Silicon$",  9, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1 );
 918   AliMedium(499, "Air gap$", 99, 0, ISXFLD, SXMGMX, 10., .1, .1, .1,   .1 );
 919 }
 920
 921 //_____________________________________________________________________________
 922 struct Bin {
 923    const AliTPCdigit *dig;
 924    int idx;
 925    Bin() {dig=0; idx=-1;}
 926 };
 927
 928 struct PreCluster : public AliTPCcluster {
 929    const AliTPCdigit* summit;
 930    int idx;
 931    int cut;
 932    int npeaks;
 933    PreCluster() : AliTPCcluster() {cut=npeaks=0;}
 934 };
 935
 936 //_____________________________________________________________________________
 937 static void FindCluster(int i, int j, Bin bins[][MAXTPCTBK+2], PreCluster &c)
 938 {
 939   //
 940   // Find clusters
 941   //
 942   Bin& b=bins[i][j];
 943   int q=b.dig->fSignal;
 944
 945   if (q<0) { // digit is at the edge of the pad row
 946     q=-q;
 947     c.cut=1;
 948   }
 949   if (b.idx >= 0 && b.idx != c.idx) {
 950     c.idx=b.idx;
 951     c.npeaks++;
 952   }
 953
 954   if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
 955
 956   c.fY += i*q;
 957   c.fZ += j*q;
 958   c.fSigmaY2 += i*i*q;
 959   c.fSigmaZ2 += j*j*q;
 960   c.fQ += q;
 961
 962   b.dig = 0;  b.idx = c.idx;
 963
 964   if (bins[i-1][j].dig) FindCluster(i-1,j,bins,c);
 965   if (bins[i][j-1].dig) FindCluster(i,j-1,bins,c);
 966   if (bins[i+1][j].dig) FindCluster(i+1,j,bins,c);
 967   if (bins[i][j+1].dig) FindCluster(i,j+1,bins,c);
 968
 969 }
 970
 971 //_____________________________________________________________________________
 972 void AliTPC::Digits2Clusters()
 973 {
 974   //
 975   // simple TPC cluster finder from digits.
 976   //
 977   const Int_t MAX_PAD=200+2, MAX_BUCKET=MAXTPCTBK+2;
 978   const Int_t Q_min=200;//75;
 979
 980   TTree *t=gAlice->TreeD();
 981   t->GetBranch("TPC")->SetAddress(&fDigits);
 982   Int_t sectors_by_rows=(Int_t)t->GetEntries();
 983
 984   int ncls=0;
 985
 986   for (Int_t n=0; n<sectors_by_rows; n++) {
 987     if (!t->GetEvent(n)) continue;
 988     Bin bins[MAX_PAD][MAX_BUCKET];
 989     AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0);
 990     Int_t nsec=dig->fSector, nrow=dig->fPadRow;
 991     Int_t ndigits=fDigits->GetEntriesFast();
 992
 993     int npads;  int sign_z;
 994     if (nsec<25) {
 995       sign_z=(nsec<13) ? 1 : -1;
 996       npads=npads_low[nrow-1];
 997     } else {
 998       sign_z=(nsec<49) ? 1 : -1;
 999       npads=npads_up[nrow-1];
1000     }
1001
1002     int ndig;
1003     for (ndig=0; ndig<ndigits; ndig++) {
1004       dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
1005       int i=dig->fPad, j=dig->fTime;
1006       if (dig->fSignal >= THRESHOLD) bins[i][j].dig=dig;
1007       if (i==1 || i==npads || j==1 || j==MAXTPCTBK) dig->fSignal*=-1;
1008     }
1009
1010     int ncl=0;
1011     int i,j;
1012     for (i=1; i<MAX_PAD-1; i++) {
1013       for (j=1; j<MAX_BUCKET-1; j++) {
1014         if (bins[i][j].dig == 0) continue;
1015         PreCluster c; c.summit=bins[i][j].dig; c.idx=ncls;
1016         FindCluster(i, j, bins, c);
1017         //if (c.fQ <= Q_min) continue; //noise cluster
1018         c.fY /= c.fQ;
1019         c.fZ /= c.fQ;
1020         c.fSigmaY2 = c.fSigmaY2/c.fQ - c.fY*c.fY + 1./12.;
1021         c.fSigmaZ2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ + 1./12.;
1022         c.fSigmaY2 *= pad_pitch_w*pad_pitch_w;
1023         c.fSigmaZ2 *= z_end/MAXTPCTBK*z_end/MAXTPCTBK;
1024         c.fSigmaY2 *= 0.022*8;
1025         c.fSigmaZ2 *= 0.068*4;
1026         c.fY = (c.fY - 0.5 - 0.5*npads)*pad_pitch_w;
1027         c.fZ = z_end/MAXTPCTBK*c.fZ;
1028         c.fZ -= 3.*fwhm/2.35482*v_drift; // PASA delay
1029         c.fZ = sign_z*(z_end - c.fZ);
1030         c.fSector=nsec-1;
1031         c.fPadRow=nrow-1;
1032         c.fTracks[0]=c.summit->fTracks[0];
1033         c.fTracks[1]=c.summit->fTracks[1];
1034         c.fTracks[2]=c.summit->fTracks[2];
1035
1036         if (c.cut) {
1037           c.fSigmaY2 *= 25.;
1038           c.fSigmaZ2 *= 4.;
1039         }
1040
1041         AddCluster(c); ncls++; ncl++;
1042       }
1043     }
1044
1045     for (ndig=0; ndig<ndigits; ndig++) {
1046       dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
1047       if (TMath::Abs(dig->fSignal) >= THRESHOLD/3)
1048         bins[dig->fPad][dig->fTime].dig=dig;
1049     }
1050
1051     for (i=1; i<MAX_PAD-1; i++) {
1052       for (j=1; j<MAX_BUCKET-1; j++) {
1053         if (bins[i][j].dig == 0) continue;
1054         PreCluster c; c.summit=bins[i][j].dig; c.idx=ncls;
1055         FindCluster(i, j, bins, c);
1056         if (c.fQ <= Q_min) continue; //noise cluster
1057         if (c.npeaks>1) continue;    //overlapped cluster
1058         c.fY /= c.fQ;
1059         c.fZ /= c.fQ;
1060         c.fSigmaY2 = c.fSigmaY2/c.fQ - c.fY*c.fY + 1./12.;
1061         c.fSigmaZ2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ + 1./12.;
1062         c.fSigmaY2 *= pad_pitch_w*pad_pitch_w;
1063         c.fSigmaZ2 *= z_end/MAXTPCTBK*z_end/MAXTPCTBK;
1064         c.fSigmaY2 *= 0.022*4;
1065         c.fSigmaZ2 *= 0.068*4;
1066         c.fY = (c.fY - 0.5 - 0.5*npads)*pad_pitch_w;
1067         c.fZ = z_end/MAXTPCTBK*c.fZ;
1068         c.fZ -= 3.*fwhm/2.35482*v_drift; // PASA delay
1069         c.fZ = sign_z*(z_end - c.fZ);
1070         c.fSector=nsec-1;
1071         c.fPadRow=nrow-1;
1072         c.fTracks[0]=c.summit->fTracks[0];
1073         c.fTracks[1]=c.summit->fTracks[1];
1074         c.fTracks[2]=c.summit->fTracks[2];
1075
1076         if (c.cut) {
1077           c.fSigmaY2 *= 25.;
1078           c.fSigmaZ2 *= 4.;
1079         }
1080
1081         if (c.npeaks==0) {AddCluster(c); ncls++; ncl++;}
1082         else {
1083           new ((*fClusters)[c.idx]) AliTPCcluster(c);
1084         }
1085       }
1086     }
1087
1088     cerr<<"sector, row, digits, clusters: "
1089         <<nsec<<' '<<nrow<<' '<<ndigits<<' '<<ncl<<"                  \r";
1090
1091     fDigits->Clear();
1092
1093   }
1094 }
1095
1096 //_____________________________________________________________________________
1097 void AliTPC::ElDiff(Float_t *xyz)
1098 {
1099   //
1100   // calculates the diffusion of a single electron
1101   //
1102
1103   Float_t driftl;
1104   //
1105   Float_t z0=xyz[2];
1106   driftl=z_end-TMath::Abs(xyz[2]);
1107   if(driftl<0.01) driftl=0.01;
1108   driftl=TMath::Sqrt(driftl);
1109   Float_t sig_t = driftl*diff_t;
1110   Float_t sig_l = driftl*diff_l;
1111   //
1112   xyz[0]=gRandom->Gaus(xyz[0],sig_t);
1113   xyz[1]=gRandom->Gaus(xyz[1],sig_t);
1114   xyz[2]=gRandom->Gaus(xyz[2],sig_l);
1115   //
1116   if (TMath::Abs(xyz[2])>z_end){
1117     xyz[2]=z_end*TMath::Sign(1.,(double) z0);
1118   }
1119   if(xyz[2]*z0 < 0.){
1120     xyz[2]=0.0001*TMath::Sign(1.,(double) z0);
1121   }
1122 }
1123
1124 //_____________________________________________________________________________
1125 void AliTPC::Hits2Clusters()
1126 {
1127   //
1128   // TPC simple cluster generator from hits
1129   // obtained from the TPC Fast Simulator
1130   //
1131   Float_t sigma_rphi,sigma_z,cl_rphi,cl_z;
1132   //
1133   GParticle *particle; // pointer to a given particle
1134   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1135   TClonesArray *Particles; //pointer to the particle list
1136   Int_t sector,nhits;
1137   Int_t ipart;
1138   Float_t xyz[5];
1139   Float_t pl,pt,tanth,rpad,ratio;
1140   Float_t rot_angle;
1141   Float_t cph,sph;
1142
1143   //---------------------------------------------------------------
1144   //  Get the access to the tracks
1145   //---------------------------------------------------------------
1146
1147   TTree *TH = gAlice->TreeH();
1148   Stat_t ntracks = TH->GetEntries();
1149
1150   //------------------------------------------------------------
1151   // Loop over all sectors (72 sectors)
1152   // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0
1153   // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0
1154   //
1155   // First cluster for sector 1 starts at "0"
1156   //------------------------------------------------------------
1157
1158
1159   fClustersIndex[0] = 0;
1160
1161   //
1162   for(Int_t isec=1;isec<fNsectors+1;isec++){
1163     //
1164     if(isec < 25){
1165       rot_angle = (isec < 13) ? (isec-1)*alpha_low : (isec-13)*alpha_low;
1166     }
1167     else {
1168       rot_angle = (isec < 49) ? (isec-25)*alpha_up : (isec-49)*alpha_up;
1169     }
1170
1171     cph=TMath::Cos(rot_angle);
1172     sph=TMath::Sin(rot_angle);
1173
1174     //------------------------------------------------------------
1175     // Loop over tracks
1176     //------------------------------------------------------------
1177
1178     for(Int_t track=0;track<ntracks;track++){
1179       ResetHits();
1180       TH->GetEvent(track);
1181       //
1182       //  Get number of the TPC hits and a pointer
1183       //  to the particles
1184       //
1185       nhits=fHits->GetEntriesFast();
1186       Particles=gAlice->Particles();
1187       //
1188       // Loop over hits
1189       //
1190       for(Int_t hit=0;hit<nhits;hit++){
1191         tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
1192         sector=tpcHit->fSector; // sector number
1193         if(sector != isec) continue; //terminate iteration
1194         ipart=tpcHit->fTrack;
1195         particle=(GParticle*)Particles->UncheckedAt(ipart);
1196         pl=particle->GetPz();
1197         pt=particle->GetPT();
1198         if(pt < 1.e-9) pt=1.e-9;
1199         tanth=pl/pt;
1200         tanth = TMath::Abs(tanth);
1201         rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY);
1202         ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
1203
1204         //   space-point resolutions
1205
1206         sigma_rphi=SigmaY2(rpad,tanth,pt);
1207         sigma_z   =SigmaZ2(rpad,tanth   );
1208
1209         //   cluster widths
1210
1211         cl_rphi=ac_rphi-bc_rphi*rpad*tanth+cc_rphi*ratio*ratio;
1212         cl_z=ac_z-bc_z*rpad*tanth+cc_z*tanth*tanth;
1213
1214         // temporary protection
1215
1216         if(sigma_rphi < 0.) sigma_rphi=0.4e-3;
1217         if(sigma_z < 0.) sigma_z=0.4e-3;
1218         if(cl_rphi < 0.) cl_rphi=2.5e-3;
1219         if(cl_z < 0.) cl_z=2.5e-5;
1220
1221         //
1222
1223         //
1224         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
1225         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
1226         //
1227         Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph;
1228         Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph;
1229         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigma_rphi));   // y
1230         xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigma_z)); // z
1231         xyz[2]=tpcHit->fQ;                                     // q
1232         xyz[3]=sigma_rphi;                                     // fSigmaY2
1233         xyz[4]=sigma_z;                                        // fSigmaZ2
1234
1235         //find row number
1236         int row;
1237         if (xprim > 0.5*(pad_row_up[0]+pad_row_low[nrow_low-1])) {
1238           for (row=0; row<nrow_up; row++) if (xprim < pad_row_up[row]) break;
1239         } else {
1240           for (row=0; row<nrow_low; row++) if (xprim < pad_row_low[row]) break;
1241         }
1242
1243         // and finally add the cluster
1244         Int_t tracks[5]={tpcHit->fTrack+1, 0, 0, sector-1, row-1};
1245         AddCluster(xyz,tracks);
1246
1247       } // end of loop over hits
1248     }   // end of loop over tracks
1249
1250     fClustersIndex[isec] = fNclusters; // update clusters index
1251
1252   } // end of loop over sectors
1253
1254   fClustersIndex[fNsectors]--; // set end of the clusters buffer
1255
1256 }
1257
1258 //_____________________________________________________________________________
1259 void AliTPC::Hits2Digits()
1260 {
1261   //
1262   // TPC conversion from hits to digits.
1263   //
1264
1265   Int_t i;
1266   //
1267   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1268   //
1269   Float_t xyz[3];
1270   Float_t rot_angle;
1271   //-------------------------------------------------------
1272   //  Get the access to the tracks
1273   //-------------------------------------------------------
1274   TTree *TH = gAlice->TreeH();
1275   Stat_t ntracks = TH->GetEntries();
1276
1277   //----------------------------------------------------
1278   // Loop over all sectors (72 sectors)
1279   // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0
1280   // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0
1281   //----------------------------------------------------
1282
1283   for(Int_t isec=1;isec<fNsectors+1;isec++){
1284
1285     //
1286     printf("*** Processing sector number %d ***\n",isec);
1287
1288     if(isec < 25){
1289       rot_angle = (isec < 13) ? (isec-1)*alpha_low : (isec-13)*alpha_low;
1290     }
1291     else {
1292       rot_angle = (isec < 49) ? (isec-25)*alpha_up : (isec-49)*alpha_up;
1293     }
1294
1295     Int_t nrows = (isec<25) ? 23 : 52;
1296
1297
1298     Float_t cph=TMath::Cos(rot_angle);
1299     Float_t sph=TMath::Sin(rot_angle);
1300
1301
1302
1303     //----------------------------------------------
1304     // Create TObjArray-s, one for each row
1305     //----------------------------------------------
1306
1307     TObjArray **row = new TObjArray* [nrows];
1308     for(i=0; i<nrows; i++){
1309       row[i] = new TObjArray;
1310     }
1311
1312     //----------------------------------------------
1313     // Loop over tracks
1314     //----------------------------------------------
1315     for(Int_t track=0;track<ntracks;track++){
1316       ResetHits();
1317       TH->GetEvent(track);
1318
1319       //------------------------------------------------
1320       //  Get number of the TPC hits and a pointer
1321       //  to the particles
1322       //------------------------------------------------
1323       Int_t nhits=fHits->GetEntriesFast();
1324       if(nhits == 0) continue;
1325       //-----------------------------------------------
1326       //  Create vectors for storing the track information,
1327       //  one vector per track per row,
1328       //  first element is a track number and then
1329       //  there are (y,z) pairs * number of electrons
1330       //----------------------------------------------
1331
1332       TVector **tr = new TVector* [nrows];
1333       Int_t *n_of_electrons= new int [nrows]; // electron counter
1334       for(i=0;i<nrows;i++){
1335         tr[i] = new TVector(241); // 120 electrons initialy
1336         n_of_electrons[i]=0;
1337       }
1338       //-----------------------------------------------------
1339       // Loop over hits
1340       //------------------------------------------------------
1341       for(Int_t hit=0;hit<nhits;hit++){
1342         tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
1343         Int_t sector=tpcHit->fSector; // sector number
1344         if(sector != isec) continue; //terminate iteration
1345
1346         xyz[0]=tpcHit->fX;
1347         xyz[1]=tpcHit->fY;
1348         xyz[2]=tpcHit->fZ;
1349         Int_t QI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
1350
1351         //-----------------------------------------------
1352         //  Rotate the electron cluster to sector 1,13,25,49
1353         //-----------------------------------------------
1354         Float_t xprim=xyz[0]*cph+xyz[1]*sph;
1355         Float_t yprim=-xyz[0]*sph+xyz[1]*cph;
1356         Float_t zprim=xyz[2];
1357
1358         //-------------------------------------
1359         // Loop over electrons
1360         //-------------------------------------
1361         for(Int_t nel=0;nel<QI;nel++){
1362           xyz[0]=xprim;  //
1363           xyz[1]=yprim;  // Keep the initial cluster position!
1364           xyz[2]=zprim;  //
1365
1366           ElDiff(xyz); // Appply the diffusion
1367
1368           Float_t row_first;
1369           Int_t row_number;
1370           row_first = (isec<25) ? pad_row_low[0] : pad_row_up[0];
1371
1372           row_number=(Int_t) ((xyz[0]-row_first+0.5*pad_pitch_l)/pad_pitch_l);
1373
1374           // Check if out of range
1375
1376           if((xyz[0]-row_first+0.5*pad_pitch_l) < 0
1377              || row_number > (nrows-1)) continue;
1378
1379           n_of_electrons[row_number]++;
1380
1381           //
1382           // Expand vector if necessary
1383           //
1384           if(n_of_electrons[row_number]>120){
1385             Int_t range = tr[row_number]->GetNrows();
1386             if(n_of_electrons[row_number] > (range-1)/2){
1387               tr[row_number]->ResizeTo(range+30); // Add 15 electrons
1388             }
1389           }
1390
1391           //---------------------------------
1392           //  E x B effect at the wires
1393           //---------------------------------
1394           Int_t nw;
1395           nw=(nwires+1)/2;
1396           Float_t xx,dx;
1397           for (Int_t nwire=1;nwire<=nwires;nwire++){
1398             xx=(nwire-nw)*ww_pitch+
1399               ((isec<13) ? pad_row_low[row_number]:pad_row_up[row_number]);
1400             dx=xx-xyz[0];
1401             if(TMath::Abs(dx) < 0.5*ww_pitch) {
1402               xyz[1]=dx*omega_tau+xyz[1];
1403               break;
1404             }
1405           } // end of loop over the wires
1406
1407           TVector &v = *tr[row_number];
1408           Int_t idx = 2*n_of_electrons[row_number]-1;
1409           v(idx)=xyz[1];
1410           v(idx+1)=xyz[2];
1411         } // end of loop over electrons
1412       } // end of loop over hits
1413
1414       //
1415       //  The track is finished
1416       //
1417       int trk=((AliTPChit*)fHits->UncheckedAt(0))->fTrack; //Y.Belikov
1418       for(i=0;i<nrows;i++){
1419         TVector &v = *tr[i];
1420         if(n_of_electrons[i] >0) {
1421           //        v(0)=(Float_t)(track+1); // track number
1422           v(0)=(Float_t)(trk+1); // Y.Belikov
1423           tr[i]->ResizeTo(2*n_of_electrons[i]+1); // shrink if necessary
1424           row[i]->Add(tr[i]); // add to the row-array
1425         }
1426         else{
1427           delete tr[i]; // delete TVector if empty
1428         }
1429       }
1430       delete [] tr; // delete track pointers
1431       delete [] n_of_electrons; // delete n_of_electrons array
1432     } //end of loop over tracks
1433     //---------------------------------------------------
1434     //  Digitize the sector data, row by row
1435     //---------------------------------------------------
1436
1437     printf("*** Digitizing sector number %d ***\n",isec);
1438
1439     for(i=0;i<nrows;i++){
1440
1441       DigSignal(isec,i,row[i]);
1442       gAlice->TreeD()->Fill();
1443       int ndig=fDigits->GetEntriesFast();
1444       printf("*** Sector, row, digits %d %d %d ***\n",isec,i+1,ndig);
1445
1446       ResetDigits();
1447
1448       row[i]->Delete(); // delete objects in this array
1449       delete row[i]; // delete the TObjArray itselves
1450
1451     } // end of digitization
1452
1453     delete [] row; // delete vectors of pointers to sectors
1454
1455   } //end of loop over sectors
1456
1457 }
1458
1459 //_____________________________________________________________________________
1460 void AliTPC::Init()
1461 {
1462   //
1463   // Initialise TPC detector after definition of geometry
1464   //
1465   Int_t i;
1466   //
1467   printf("\n");
1468   for(i=0;i<35;i++) printf("*");
1469   printf(" TPC_INIT ");
1470   for(i=0;i<35;i++) printf("*");
1471   printf("\n");
1472   //
1473   for(i=0;i<80;i++) printf("*");
1474   printf("\n");
1475 }
1476
1477 //_____________________________________________________________________________
1478 void AliTPC::MakeBranch(Option_t* option)
1479 {
1480   //
1481   // Create Tree branches for the TPC.
1482   //
1483   Int_t buffersize = 4000;
1484   char branchname[10];
1485   sprintf(branchname,"%s",GetName());
1486
1487   AliDetector::MakeBranch(option);
1488
1489   char *D = strstr(option,"D");
1490
1491   if (fDigits   && gAlice->TreeD() && D) {
1492     gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
1493     printf("Making Branch %s for digits\n",branchname);
1494   }
1495
1496   char *R = strstr(option,"R");
1497
1498   if (fClusters && gAlice->TreeR() && R) {
1499     gAlice->TreeR()->Branch(branchname,&fClusters, buffersize);
1500     printf("Making Branch %s for Clusters\n",branchname);
1501   }
1502 }
1503
1504 //_____________________________________________________________________________
1505 void AliTPC::ResetDigits()
1506 {
1507   //
1508   // Reset number of digits and the digits array for this detector
1509   // reset clusters
1510   //
1511   fNdigits   = 0;
1512   if (fDigits)   fDigits->Clear();
1513   fNclusters = 0;
1514   if (fClusters) fClusters->Clear();
1515 }
1516
1517 //_____________________________________________________________________________
1518 void AliTPC::SetSecAL(Int_t sec)
1519 {
1520   //
1521   // Activate/deactivate selection for lower sectors
1522   //
1523   fSecAL = sec;
1524 }
1525
1526 //_____________________________________________________________________________
1527 void AliTPC::SetSecAU(Int_t sec)
1528 {
1529   //
1530   // Activate/deactivate selection for upper sectors
1531   //
1532   fSecAU = sec;
1533 }
1534
1535 //_____________________________________________________________________________
1536 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
1537 {
1538   //
1539   // Select active lower sectors
1540   //
1541   fSecLows[0] = s1;
1542   fSecLows[1] = s2;
1543   fSecLows[2] = s3;
1544   fSecLows[3] = s4;
1545   fSecLows[4] = s5;
1546   fSecLows[5] = s6;
1547 }
1548
1549 //_____________________________________________________________________________
1550 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
1551                        Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
1552                        Int_t s11 , Int_t s12)
1553 {
1554   //
1555   // Select active upper sectors
1556   //
1557   fSecUps[0] = s1;
1558   fSecUps[1] = s2;
1559   fSecUps[2] = s3;
1560   fSecUps[3] = s4;
1561   fSecUps[4] = s5;
1562   fSecUps[5] = s6;
1563   fSecUps[6] = s7;
1564   fSecUps[7] = s8;
1565   fSecUps[8] = s9;
1566   fSecUps[9] = s10;
1567   fSecUps[10] = s11;
1568   fSecUps[11] = s12;
1569 }
1570
1571 //_____________________________________________________________________________
1572 void AliTPC::SetSens(Int_t sens)
1573 {
1574   fSens = sens;
1575 }
1576
1577 //_____________________________________________________________________________
1578 void AliTPC::Streamer(TBuffer &R__b)
1579 {
1580   //
1581   // Stream an object of class AliTPC.
1582   //
1583    if (R__b.IsReading()) {
1584       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1585       AliDetector::Streamer(R__b);
1586       if (R__v < 2) return;
1587       R__b >> fNsectors;
1588       R__b >> fNclusters;
1589       R__b >> fNtracks;
1590       fClustersIndex = new Int_t[fNsectors+1];
1591       fDigitsIndex   = new Int_t[fNsectors+1];
1592    } else {
1593       R__b.WriteVersion(AliTPC::IsA());
1594       AliDetector::Streamer(R__b);
1595       R__b << fNsectors;
1596       R__b << fNclusters;
1597       R__b << fNtracks;
1598    }
1599 }
1600
1601 ClassImp(AliTPCcluster)
1602
1603 //_____________________________________________________________________________
1604 AliTPCcluster::AliTPCcluster(Float_t *hits, Int_t *lab)
1605 {
1606   //
1607   // Creates a simulated cluster for the TPC
1608   //
1609   fTracks[0]  = lab[0];
1610   fTracks[1]  = lab[1];
1611   fTracks[2]  = lab[2];
1612   fSector     = lab[3];
1613   fPadRow     = lab[4];
1614   fY          = hits[0];
1615   fZ          = hits[1];
1616   fQ          = hits[2];
1617   fSigmaY2    = hits[3];
1618   fSigmaZ2    = hits[4];
1619 }
1620
1621 //_____________________________________________________________________________
1622 void AliTPCcluster::GetXYZ(Double_t& x, Double_t& y, Double_t& z) const
1623 {
1624   //
1625   // Returns centroid for of a simulated TPC cluster
1626   //
1627   Double_t alpha,xrow;
1628   if (fSector<24) {
1629     alpha=(fSector<12) ? fSector*alpha_low : (fSector-12)*alpha_low;
1630     xrow=pad_row_low[fPadRow];
1631   } else {
1632     alpha=(fSector<48) ? (fSector-24)*alpha_up : (fSector-48)*alpha_up;
1633     xrow=pad_row_up[fPadRow];
1634   }
1635   x=xrow*cos(alpha) - fY*sin(alpha);
1636   y=xrow*sin(alpha) + fY*cos(alpha);
1637   z=fZ;
1638 }
1639
1640 ClassImp(AliTPCdigit)
1641
1642 //_____________________________________________________________________________
1643 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
1644   AliDigit(tracks)
1645 {
1646   //
1647   // Creates a TPC digit object
1648   //
1649   fSector     = digits[0];
1650   fPadRow     = digits[1];
1651   fPad        = digits[2];
1652   fTime       = digits[3];
1653   fSignal     = digits[4];
1654 }
1655
1656
1657 ClassImp(AliTPChit)
1658
1659 //_____________________________________________________________________________
1660 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
1661 AliHit(shunt,track)
1662 {
1663   //
1664   // Creates a TPC hit object
1665   //
1666   fSector     = vol[0];
1667   fPadRow     = vol[1];
1668   fX          = hits[0];
1669   fY          = hits[1];
1670   fZ          = hits[2];
1671   fQ          = hits[3];
1672 }
1673
1674
1675 ClassImp(AliTPCtrack)
1676
1677 //_____________________________________________________________________________
1678 AliTPCtrack::AliTPCtrack(Float_t *hits)
1679 {
1680   //
1681   // Default creator for a TPC reconstructed track object
1682   //
1683   ref=hits[0]; // This is dummy code !
1684 }
1685
1686 AliTPCtrack::AliTPCtrack(const AliTPCcluster& c,const TVector& xx,
1687                          const TMatrix& CC):
1688   x(xx),C(CC),clusters(MAX_CLUSTER)
1689 {
1690   //
1691   // Standard creator for a TPC reconstructed track object
1692   //
1693   chi2=0.;
1694   int sec=c.fSector, row=c.fPadRow;
1695   if (sec<24) {
1696     fAlpha=(sec%12)*alpha_low; ref=pad_row_low[0] + row*pad_pitch_l;
1697   } else {
1698     fAlpha=(sec%24)*alpha_up;  ref=pad_row_up[0]  + row*pad_pitch_l;
1699   }
1700   clusters.AddLast((AliTPCcluster*)(&c));
1701 }
1702
1703 //_____________________________________________________________________________
1704 AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : x(t.x), C(t.C),
1705   clusters(t.clusters.GetEntriesFast())
1706 {
1707   //
1708   // Copy creator for a TPC reconstructed track
1709   //
1710   ref=t.ref;
1711   chi2=t.chi2;
1712   fAlpha=t.fAlpha;
1713   int n=t.clusters.GetEntriesFast();
1714   for (int i=0; i<n; i++) clusters.AddLast(t.clusters.UncheckedAt(i));
1715 }
1716
1717 //_____________________________________________________________________________
1718 Double_t AliTPCtrack::GetY(Double_t xk) const
1719 {
1720   //
1721   //
1722   //
1723   Double_t c2=x(2)*xk - x(3);
1724   if (TMath::Abs(c2) >= 0.999) {
1725     if (*this>10) cerr<<*this<<" AliTPCtrack warning: No y for this x !\n";
1726     return 0.;
1727   }
1728   Double_t c1=x(2)*ref - x(3);
1729   Double_t r1=sqrt(1.-c1*c1), r2=sqrt(1.-c2*c2);
1730   Double_t dx=xk-ref;
1731   return x(0) + dx*(c1+c2)/(r1+r2);
1732 }
1733
1734 //_____________________________________________________________________________
1735 int AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
1736 {
1737   //
1738   // Propagate a TPC reconstructed track
1739   //
1740   if (TMath::Abs(x(2)*xk - x(3)) >= 0.999) {
1741     if (*this>3) cerr<<*this<<" AliTPCtrack warning: Propagation failed !\n";
1742     return 0;
1743   }
1744
1745   Double_t x1=ref, x2=x1+0.5*(xk-x1), dx=x2-x1, y1=x(0), z1=x(1);
1746   Double_t c1=x(2)*x1 - x(3), r1=sqrt(1.- c1*c1);
1747   Double_t c2=x(2)*x2 - x(3), r2=sqrt(1.- c2*c2);
1748
1749   x(0) += dx*(c1+c2)/(r1+r2);
1750   x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
1751
1752   TMatrix F(5,5); F.UnitMatrix();
1753   Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2;
1754   F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
1755   F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
1756   Double_t cr=c1*r2+c2*r1;
1757   F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
1758   F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
1759   F(1,4)= dx*cc/cr;
1760   TMatrix tmp(F,TMatrix::kMult,C);
1761   C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
1762
1763   ref=x2;
1764
1765   //Multiple scattering******************
1766   Double_t ey=x(2)*ref - x(3);
1767   Double_t ex=sqrt(1-ey*ey);
1768   Double_t ez=x(4);
1769   TMatrix Q(5,5); Q=0.;
1770   Q(2,2)=ez*ez+ey*ey;   Q(2,3)=-ex*ey;       Q(2,4)=-ex*ez;
1771   Q(3,2)=Q(2,3);        Q(3,3)= ez*ez+ex*ex; Q(3,4)=-ey*ez;
1772   Q(4,2)=Q(2,4);        Q(4,3)= Q(3,4);      Q(4,4)=1.;
1773
1774   F=0;
1775   F(2,2)=-x(2)*ex;          F(2,3)=-x(2)*ey;
1776   F(3,2)=-ex*(x(2)*ref-ey); F(3,3)=-(1.+ x(2)*ref*ey - ey*ey);
1777   F(4,2)=-ez*ex;            F(4,3)=-ez*ey;           F(4,4)=1.;
1778
1779   tmp.Mult(F,Q);
1780   Q.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
1781
1782   Double_t p2=GetPt()*GetPt()*(1.+x(4)*x(4));
1783   Double_t beta2=p2/(p2 + pm*pm);
1784   Double_t d=sqrt((x1-ref)*(x1-ref)+(y1-x(0))*(y1-x(0))+(z1-x(1))*(z1-x(1)));
1785   d*=2.;
1786   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
1787   Q*=theta2;
1788   C+=Q;
1789
1790   //Energy losses************************
1791   Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho;
1792   if (x1 < x2) dE=-dE;
1793   x(2)*=(1.- sqrt(p2+pm*pm)/p2*dE);
1794
1795   x1=ref; x2=xk; y1=x(0); z1=x(1);
1796   c1=x(2)*x1 - x(3); r1=sqrt(1.- c1*c1);
1797   c2=x(2)*x2 - x(3); r2=sqrt(1.- c2*c2);
1798
1799   x(0) += dx*(c1+c2)/(r1+r2);
1800   x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
1801
1802   F.UnitMatrix();
1803   rr=r1+r2; cc=c1+c2; xx=x1+x2;
1804   F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
1805   F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
1806   cr=c1*r2+c2*r1;
1807   F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
1808   F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
1809   F(1,4)= dx*cc/cr;
1810   tmp.Mult(F,C);
1811   C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
1812
1813   ref=x2;
1814
1815   return 1;
1816 }
1817
1818 //_____________________________________________________________________________
1819 void AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm)
1820 {
1821   //
1822   // Propagate a reconstructed track from the vertex
1823   //
1824   Double_t c=x(2)*ref - x(3);
1825   Double_t tgf=-x(3)/(x(2)*x(0) + sqrt(1-c*c));
1826   Double_t snf=tgf/sqrt(1.+ tgf*tgf);
1827   Double_t xv=(x(3)+snf)/x(2);
1828   PropagateTo(xv,x0,rho,pm);
1829 }
1830
1831 //_____________________________________________________________________________
1832 void AliTPCtrack::Update(const AliTPCcluster *c, Double_t chisq)
1833 {
1834   //
1835   // Update statistics for a reconstructed TPC track
1836   //
1837   TMatrix H(2,5); H.UnitMatrix();
1838   TMatrix Ht(TMatrix::kTransposed,H);
1839   TVector m(2);   m(0)=c->fY; m(1)=c->fZ;
1840   TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
1841
1842   TMatrix tmp(H,TMatrix::kMult,C);
1843   TMatrix R(tmp,TMatrix::kMult,Ht); R+=V;
1844
1845   Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
1846   R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
1847   R(1,0)*=-1; R(0,1)=R(1,0);
1848   R*=1./det;
1849
1850   //R.Invert();
1851
1852   TMatrix K(C,TMatrix::kMult,Ht); K*=R;
1853
1854   TVector savex=x;
1855   x*=H; x-=m; x*=-1; x*=K; x+=savex;
1856   if (TMath::Abs(x(2)*ref-x(3)) >= 0.999) {
1857     if (*this>3) cerr<<*this<<" AliTPCtrack warning: Filtering failed !\n";
1858     x=savex;
1859     return;
1860   }
1861
1862   TMatrix saveC=C;
1863   C.Mult(K,tmp); C-=saveC; C*=-1;
1864
1865   clusters.AddLast((AliTPCcluster*)c);
1866   chi2 += chisq;
1867 }
1868
1869 //_____________________________________________________________________________
1870 int AliTPCtrack::Rotate(Double_t alpha)
1871 {
1872   //
1873   // Rotate a reconstructed TPC track
1874   //
1875   fAlpha += alpha;
1876
1877   Double_t x1=ref, y1=x(0);
1878   Double_t ca=cos(alpha), sa=sin(alpha);
1879   Double_t r1=x(2)*ref - x(3);
1880
1881   ref = x1*ca + y1*sa;
1882   x(0)=-x1*sa + y1*ca;
1883   x(3)=x(3)*ca + (x(2)*y1 + sqrt(1.- r1*r1))*sa;
1884
1885   Double_t r2=x(2)*ref - x(3);
1886   if (TMath::Abs(r2) >= 0.999) {
1887     if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !\n";
1888     return 0;
1889   }
1890
1891   Double_t y0=x(0) + sqrt(1.- r2*r2)/x(2);
1892   if ((x(0)-y0)*x(2) >= 0.) {
1893     if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !!!\n";
1894     return 0;
1895   }
1896
1897   TMatrix F(5,5); F.UnitMatrix();
1898   F(0,0)=ca;
1899   F(3,0)=x(2)*sa;
1900   F(3,2)=(y1 - r1*x1/sqrt(1.- r1*r1))*sa;
1901   F(3,3)= ca + sa*r1/sqrt(1.- r1*r1);
1902   TMatrix tmp(F,TMatrix::kMult,C);
1903   // Double_t dy2=C(0,0);
1904   C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
1905   // C(0,0)+=dy2*sa*sa*r1*r1/(1.- r1*r1);
1906   // C(1,1)+=dy2*sa*sa*x(4)*x(4)/(1.- r1*r1);
1907
1908   return 1;
1909 }
1910
1911 //_____________________________________________________________________________
1912 void AliTPCtrack::UseClusters() const
1913 {
1914   //
1915   //
1916   //
1917   int num_of_clusters=clusters.GetEntriesFast();
1918   for (int i=0; i<num_of_clusters; i++) {
1919     //if (i<=14) continue;
1920     AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i);
1921     c->Use();
1922   }
1923 }
1924
1925 //_____________________________________________________________________________
1926 Double_t AliTPCtrack::GetPredictedChi2(const AliTPCcluster *c) const
1927 {
1928   //
1929   // Calculate chi2 for a reconstructed TPC track
1930   //
1931   TMatrix H(2,5); H.UnitMatrix();
1932   TVector m(2);   m(0)=c->fY; m(1)=c->fZ;
1933   TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
1934   TVector res=x;  res*=H; res-=m; //res*=-1;
1935   TMatrix tmp(H,TMatrix::kMult,C);
1936   TMatrix R(tmp,TMatrix::kMult,TMatrix(TMatrix::kTransposed,H)); R+=V;
1937
1938   Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
1939   if (TMath::Abs(det) < 1.e-10) {
1940     if (*this>3) cerr<<*this<<" AliTPCtrack warning: Singular matrix !\n";
1941     return 1e10;
1942   }
1943   R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
1944   R(1,0)*=-1; R(0,1)=R(1,0);
1945   R*=1./det;
1946
1947   //R.Invert();
1948
1949   TVector r=res;
1950   res*=R;
1951   return r*res;
1952 }
1953
1954 //_____________________________________________________________________________
1955 int AliTPCtrack::GetLab() const
1956 {
1957   //
1958   //
1959   //
1960   int lab = 0;
1961   struct {
1962     int lab;
1963     int max;
1964   } s[MAX_CLUSTER]={{0,0}};
1965
1966   int i;
1967   int num_of_clusters=clusters.GetEntriesFast();
1968   for (i=0; i<num_of_clusters; i++) {
1969     AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i);
1970     lab=c->fTracks[0]; if (lab<0) lab=-lab;
1971     int j;
1972     for (j=0; j<MAX_CLUSTER; j++)
1973       if (s[j].lab==lab || s[j].max==0) break;
1974     s[j].lab=lab;
1975     s[j].max++;
1976   }
1977
1978   int max=0;
1979   for (i=0; i<num_of_clusters; i++)
1980     if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
1981   if (lab>0) lab--;
1982
1983   if (1.-float(max)/num_of_clusters > 0.10) return -lab;
1984
1985   if (num_of_clusters < 6) return lab;
1986
1987   max=0;
1988   for (i=1; i<=6; i++) {
1989     AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(num_of_clusters-i);
1990     if (lab != TMath::Abs(c->fTracks[0])-1
1991         && lab != TMath::Abs(c->fTracks[1])-1
1992         && lab != TMath::Abs(c->fTracks[2])-1
1993         ) max++;
1994   }
1995   if (max>3) return -lab;
1996
1997   return lab;
1998 }
1999
2000 //_____________________________________________________________________________
2001 void AliTPCtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const
2002 {
2003   //
2004   // Get reconstructed TPC track momentum
2005   //
2006   Double_t pt=0.3*FIELD/TMath::Abs(x(2))/100; // GeV/c
2007   Double_t r=x(2)*ref-x(3);
2008   Double_t y0=x(0) + sqrt(1.- r*r)/x(2);
2009   px=-pt*(x(0)-y0)*x(2);    //cos(phi);
2010   py=-pt*(x(3)-ref*x(2));   //sin(phi);
2011   pz=pt*x(4);
2012   Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha);
2013   py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha);
2014   px=tmp;
2015 }
2016
2017 //_____________________________________________________________________________
2018 //
2019 //     Classes for internal tracking use
2020 //
2021
2022 //_____________________________________________________________________________
2023 void AliTPCRow::InsertCluster(const AliTPCcluster* c)
2024 {
2025   //
2026   // Insert a cluster in the list
2027   //
2028   if (num_of_clusters==MAX_CLUSTER_PER_ROW) {
2029     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
2030   }
2031   if (num_of_clusters==0) {clusters[num_of_clusters++]=c; return;}
2032   int i=Find(c->fY);
2033   memmove(clusters+i+1 ,clusters+i,(num_of_clusters-i)*sizeof(AliTPCcluster*));
2034   clusters[i]=c; num_of_clusters++;
2035 }
2036
2037 //_____________________________________________________________________________
2038 int AliTPCRow::Find(Double_t y) const
2039 {
2040   //
2041   //
2042   //
2043   if (y <= clusters[0]->fY) return 0;
2044   if (y > clusters[num_of_clusters-1]->fY) return num_of_clusters;
2045   int b=0, e=num_of_clusters-1, m=(b+e)/2;
2046   for (; b<e; m=(b+e)/2) {
2047     if (y > clusters[m]->fY) b=m+1;
2048     else e=m;
2049   }
2050   return m;
2051 }