TPC/AliTPCFast.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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  fast TPC cluster simulation                                              //
  21 //                                                                          //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TParticle.h>
  25
  26 #include "AliRunLoader.h"
  27 #include "AliRun.h"
  28 #include "AliMC.h"
  29 #include "AliTPC.h"
  30 #include "AliTPCParam.h"
  31 #include "AliTPCClustersArray.h"
  32 #include "AliTPCClustersRow.h"
  33 #include "AliTPCcluster.h"
  34 #include "AliComplexCluster.h"
  35 #include "AliTPCFast.h"
  36
  37
  38 ClassImp(AliTPCFast)
  39
  40
  41 //_____________________________________________________________________________
  42 void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
  43 {
  44   //--------------------------------------------------------
  45   // TPC simple cluster generator from hits
  46   // obtained from the TPC Fast Simulator
  47   // The point errors are taken from the parametrization
  48   //--------------------------------------------------------
  49
  50   //-----------------------------------------------------------------
  51   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
  52   //-----------------------------------------------------------------
  53   // Adopted to Marian's cluster data structure by I.Belikov, CERN,
  54   // Jouri.Belikov@cern.ch
  55   //----------------------------------------------------------------
  56
  57   /////////////////////////////////////////////////////////////////////////////
  58   //
  59   //---------------------------------------------------------------------
  60   //   ALICE TPC Cluster Parameters
  61   //--------------------------------------------------------------------
  62
  63
  64
  65   // Cluster width in rphi
  66   const Float_t kACrphi=0.18322;
  67   const Float_t kBCrphi=0.59551e-3;
  68   const Float_t kCCrphi=0.60952e-1;
  69   // Cluster width in z
  70   const Float_t kACz=0.19081;
  71   const Float_t kBCz=0.55938e-3;
  72   const Float_t kCCz=0.30428;
  73
  74
  75   AliLoader* loader = runLoader->GetLoader("TPCLoader");
  76   if (!loader) {
  77     Error("Hits2Clusters", "No TPC loader found");
  78     return;
  79   }
  80   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
  81   AliRun* aliRun = runLoader->GetAliRun();
  82   if (!aliRun) {
  83     Error("Hits2Clusters", "Couldn't get AliRun object");
  84     return;
  85   }
  86   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
  87   if (!tpc) {
  88     Error("Hits2Clusters", "Couldn't get TPC detector");
  89     return;
  90   }
  91   AliTPCParam* param = tpc->GetParam();
  92   if (!param) {
  93     Error("Hits2Clusters", "No TPC parameters available");
  94     return;
  95   }
  96
  97   //if(fDefaults == 0) SetDefaults();
  98
  99   Float_t sigmaRphi,sigmaZ,clRphi,clZ;
 100   //
 101   TParticle *particle; // pointer to a given particle
 102   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 103   Int_t sector;
 104   Int_t ipart;
 105   Float_t xyz[5];
 106   Float_t pl,pt,tanth,rpad,ratio;
 107   Float_t cph,sph;
 108
 109   //---------------------------------------------------------------
 110   //  Get the access to the tracks
 111   //---------------------------------------------------------------
 112
 113   TTree *tH = loader->TreeH();
 114   if (tH == 0x0)
 115    {
 116      Fatal("Hits2Clusters","Can not find TreeH in folder");
 117      return;
 118    }
 119   tpc->SetTreeAddress();
 120
 121   Stat_t ntracks = tH->GetEntries();
 122
 123   //Switch to the output file
 124
 125   if (loader->TreeR() == 0x0) loader->MakeTree("R");
 126
 127   cout<<"param->GetTitle() = "<<param->GetTitle()<<endl;
 128
 129   runLoader->CdGAFile();
 130   //param->Write(param->GetTitle());
 131
 132   AliTPCClustersArray carray;
 133   carray.Setup(param);
 134   carray.SetClusterType("AliTPCcluster");
 135   carray.MakeTree(loader->TreeR());
 136
 137   Int_t nclusters=0; //cluster counter
 138
 139   //------------------------------------------------------------
 140   // Loop over all sectors (72 sectors for 20 deg
 141   // segmentation for both lower and upper sectors)
 142   // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
 143   // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
 144   //
 145   // First cluster for sector 0 starts at "0"
 146   //------------------------------------------------------------
 147
 148   for(Int_t isec=0;isec<param->GetNSector();isec++){
 149     //MI change
 150     param->AdjustCosSin(isec,cph,sph);
 151
 152     //------------------------------------------------------------
 153     // Loop over tracks
 154     //------------------------------------------------------------
 155
 156     for(Int_t track=0;track<ntracks;track++){
 157       tpc->ResetHits();
 158       tH->GetEvent(track);
 159       //
 160       //  Get number of the TPC hits
 161       //
 162        tpcHit = (AliTPChit*)tpc->FirstHit(-1);
 163
 164       // Loop over hits
 165       //
 166        while(tpcHit){
 167
 168          if (tpcHit->fQ == 0.) {
 169            tpcHit = (AliTPChit*) tpc->NextHit();
 170            continue; //information about track (I.Belikov)
 171          }
 172         sector=tpcHit->fSector; // sector number
 173
 174        if(sector != isec){
 175          tpcHit = (AliTPChit*) tpc->NextHit();
 176          continue;
 177        }
 178         ipart=tpcHit->Track();
 179         particle=aliRun->GetMCApp()->Particle(ipart);
 180         pl=particle->Pz();
 181         pt=particle->Pt();
 182         if(pt < 1.e-9) pt=1.e-9;
 183         tanth=pl/pt;
 184         tanth = TMath::Abs(tanth);
 185         rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
 186         ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
 187
 188         //   space-point resolutions
 189
 190         sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
 191         sigmaZ   =AliTPCcluster::SigmaZ2(rpad,tanth   );
 192
 193         //   cluster widths
 194
 195         clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
 196         clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
 197
 198         // temporary protection
 199
 200         if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
 201         if(sigmaZ < 0.) sigmaZ=0.4e-3;
 202         if(clRphi < 0.) clRphi=2.5e-3;
 203         if(clZ < 0.) clZ=2.5e-5;
 204
 205         //
 206
 207         //
 208         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
 209         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
 210         //
 211         Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
 212         Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
 213         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
 214           Float_t alpha=(isec < param->GetNInnerSector()) ?
 215           param->GetInnerAngle() : param->GetOuterAngle();
 216           Float_t ymax=xprim*TMath::Tan(0.5*alpha);
 217           if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
 218         xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
 219           if (TMath::Abs(xyz[1])>param->GetZLength()) xyz[1]=tpcHit->Z();
 220         xyz[2]=sigmaRphi;                                     // fSigmaY2
 221         xyz[3]=sigmaZ;                                        // fSigmaZ2
 222         xyz[4]=tpcHit->fQ;                                    // q
 223
 224         AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
 225         if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
 226
 227         Int_t tracks[3]={tpcHit->Track(), -1, -1};
 228         AliTPCcluster cluster(tracks,xyz);
 229
 230         clrow->InsertCluster(&cluster); nclusters++;
 231
 232         tpcHit = (AliTPChit*)tpc->NextHit();
 233
 234
 235       } // end of loop over hits
 236
 237     }   // end of loop over tracks
 238
 239     Int_t nrows=param->GetNRow(isec);
 240     for (Int_t irow=0; irow<nrows; irow++) {
 241         AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
 242         if (!clrow) continue;
 243         carray.StoreRow(isec,irow);
 244         carray.ClearRow(isec,irow);
 245     }
 246
 247   } // end of loop over sectors
 248
 249   //  cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
 250   loader->WriteRecPoints("OVERWRITE");
 251
 252
 253 } // end of function
 254
 255 //_________________________________________________________________
 256 void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
 257                                           AliTPCClustersArray* clustersArray,
 258                                           Int_t isec) const
 259 {
 260   //--------------------------------------------------------
 261   //calculate exact cross point of track and given pad row
 262   //resulting values are expressed in "digit" coordinata
 263   //--------------------------------------------------------
 264
 265   //-----------------------------------------------------------------
 266   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
 267   //-----------------------------------------------------------------
 268   //
 269   if (clustersArray==0){
 270     return;
 271   }
 272   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 273   if (!loader) {
 274     Error("Hits2ExactClustersSector", "No TPC loader found");
 275     return;
 276   }
 277   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 278   AliRun* aliRun = runLoader->GetAliRun();
 279   if (!aliRun) {
 280     Error("Hits2ExactClustersSector", "Couldn't get AliRun object");
 281     return;
 282   }
 283   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
 284   if (!tpc) {
 285     Error("Hits2ExactClustersSector", "Couldn't get TPC detector");
 286     return;
 287   }
 288   AliTPCParam* param = tpc->GetParam();
 289   if (!param) {
 290     Error("Hits2ExactClustersSector", "No TPC parameters available");
 291     return;
 292   }
 293   //
 294   TParticle *particle; // pointer to a given particle
 295   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 296   //  Int_t sector,nhits;
 297   Int_t ipart;
 298   const Int_t kcmaxhits=30000;
 299   AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
 300   AliTPCFastVector & xxx = *xxxx;
 301   Int_t maxhits = kcmaxhits;
 302   //construct array for each padrow
 303   for (Int_t i=0; i<param->GetNRow(isec);i++)
 304     clustersArray->CreateRow(isec,i);
 305
 306   //---------------------------------------------------------------
 307   //  Get the access to the tracks
 308   //---------------------------------------------------------------
 309
 310   TTree *tH = loader->TreeH();
 311   if (tH == 0x0)
 312    {
 313      Fatal("Hits2Clusters","Can not find TreeH in folder");
 314      return;
 315    }
 316   tpc->SetTreeAddress();
 317
 318   Stat_t ntracks = tH->GetEntries();
 319   Int_t npart = aliRun->GetMCApp()->GetNtrack();
 320   //MI change
 321   TBranch * branch=0;
 322   if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
 323   else branch = tH->GetBranch("TPC");
 324
 325   //------------------------------------------------------------
 326   // Loop over tracks
 327   //------------------------------------------------------------
 328
 329   for(Int_t track=0;track<ntracks;track++){
 330     Bool_t isInSector=kTRUE;
 331     tpc->ResetHits();
 332      isInSector = tpc->TrackInVolume(isec,track);
 333     if (!isInSector) continue;
 334     //MI change
 335     branch->GetEntry(track); // get next track
 336     //
 337     //  Get number of the TPC hits and a pointer
 338     //  to the particles
 339     // Loop over hits
 340     //
 341     Int_t currentIndex=0;
 342     Int_t lastrow=-1;  //last writen row
 343
 344     //M.I. changes
 345
 346     tpcHit = (AliTPChit*)tpc->FirstHit(-1);
 347     while(tpcHit){
 348
 349       Int_t sector=tpcHit->fSector; // sector number
 350       if(sector != isec){
 351         tpcHit = (AliTPChit*) tpc->NextHit();
 352         continue;
 353       }
 354
 355       ipart=tpcHit->Track();
 356       if (ipart<npart) particle=aliRun->GetMCApp()->Particle(ipart);
 357
 358       //find row number
 359
 360       Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
 361       Int_t    index[3]={1,isec,0};
 362       Int_t    currentrow = param->GetPadRow(x,index) ;
 363       if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
 364       if (lastrow<0) lastrow=currentrow;
 365       if (currentrow==lastrow){
 366         if ( currentIndex>=maxhits){
 367           maxhits+=kcmaxhits;
 368           xxx.ResizeTo(4*maxhits);
 369         }
 370         xxx(currentIndex*4)=x[0];
 371         xxx(currentIndex*4+1)=x[1];
 372         xxx(currentIndex*4+2)=x[2];
 373         xxx(currentIndex*4+3)=tpcHit->fQ;
 374         currentIndex++;
 375       }
 376       else
 377         if (currentIndex>2){
 378           Float_t sumx=0;
 379           Float_t sumx2=0;
 380           Float_t sumx3=0;
 381           Float_t sumx4=0;
 382           Float_t sumy=0;
 383           Float_t sumxy=0;
 384           Float_t sumx2y=0;
 385           Float_t sumz=0;
 386           Float_t sumxz=0;
 387           Float_t sumx2z=0;
 388           Float_t sumq=0;
 389           for (Int_t index=0;index<currentIndex;index++){
 390             Float_t x,x2,x3,x4;
 391             x=x2=x3=x4=xxx(index*4);
 392             x2*=x;
 393             x3*=x2;
 394             x4*=x3;
 395             sumx+=x;
 396             sumx2+=x2;
 397             sumx3+=x3;
 398             sumx4+=x4;
 399             sumy+=xxx(index*4+1);
 400             sumxy+=xxx(index*4+1)*x;
 401             sumx2y+=xxx(index*4+1)*x2;
 402             sumz+=xxx(index*4+2);
 403             sumxz+=xxx(index*4+2)*x;
 404             sumx2z+=xxx(index*4+2)*x2;
 405             sumq+=xxx(index*4+3);
 406           }
 407           Float_t centralPad = (param->GetNPads(isec,lastrow)-1)/2;
 408           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
 409             sumx2*(sumx*sumx3-sumx2*sumx2);
 410
 411           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
 412             sumx2*(sumxy*sumx3-sumx2y*sumx2);
 413           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
 414             sumx2*(sumxz*sumx3-sumx2z*sumx2);
 415
 416           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
 417             sumx2*(sumx*sumx2y-sumx2*sumxy);
 418           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
 419             sumx2*(sumx*sumx2z-sumx2*sumxz);
 420
 421           if (TMath::Abs(det)<0.00001){
 422              tpcHit = (AliTPChit*)tpc->NextHit();
 423             continue;
 424           }
 425
 426           Float_t y=detay/det+centralPad;
 427           Float_t z=detaz/det;
 428           Float_t by=detby/det; //y angle
 429           Float_t bz=detbz/det; //z angle
 430           sumy/=Float_t(currentIndex);
 431           sumz/=Float_t(currentIndex);
 432
 433           AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
 434           if (row!=0) {
 435             AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
 436             //    AliComplexCluster cl;
 437             cl->fX=z;
 438             cl->fY=y;
 439             cl->fQ=sumq;
 440             cl->fSigmaX2=bz;
 441             cl->fSigmaY2=by;
 442             cl->fTracks[0]=ipart;
 443           }
 444           currentIndex=0;
 445           lastrow=currentrow;
 446         } //end of calculating cluster for given row
 447
 448
 449       tpcHit = (AliTPChit*)tpc->NextHit();
 450     } // end of loop over hits
 451   }   // end of loop over tracks
 452   //write padrows to tree
 453   for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
 454     clustersArray->StoreRow(isec,ii);
 455     clustersArray->ClearRow(isec,ii);
 456   }
 457   xxxx->Delete();
 458
 459 }
 460
 461
 462 void AliTPCFast::FindTrackHitsIntersection(AliRunLoader* runLoader,
 463                                            AliTPCClustersArray* clustersArray) const
 464 {
 465
 466   //
 467   //fill clones array with intersection of current point with the
 468   //middle of the row
 469   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 470   if (!loader) {
 471     Error("FindTrackHitsIntersection", "No TPC loader found");
 472     return;
 473   }
 474   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 475   AliRun* aliRun = runLoader->GetAliRun();
 476   if (!aliRun) {
 477     Error("FindTrackHitsIntersection", "Couldn't get AliRun object");
 478     return;
 479   }
 480   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
 481   if (!tpc) {
 482     Error("FindTrackHitsIntersection", "Couldn't get TPC detector");
 483     return;
 484   }
 485   AliTPCParam* param = tpc->GetParam();
 486   if (!param) {
 487     Error("FindTrackHitsIntersection", "No TPC parameters available");
 488     return;
 489   }
 490
 491   Int_t sector;
 492   Int_t ipart;
 493
 494   const Int_t kcmaxhits=30000;
 495   AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
 496   AliTPCFastVector & xxx = *xxxx;
 497   Int_t maxhits = kcmaxhits;
 498
 499   //
 500   AliTPChit * tpcHit=0;
 501   tpcHit = (AliTPChit*)tpc->FirstHit2(-1);
 502   Int_t currentIndex=0;
 503   Int_t lastrow=-1;  //last writen row
 504
 505   while (tpcHit){
 506     if (tpcHit==0) continue;
 507     sector=tpcHit->fSector; // sector number
 508     ipart=tpcHit->Track();
 509
 510     //find row number
 511
 512     Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
 513     Int_t    index[3]={1,sector,0};
 514     Int_t    currentrow = param->GetPadRow(x,index) ;
 515     if (currentrow<0) continue;
 516     if (lastrow<0) lastrow=currentrow;
 517     if (currentrow==lastrow){
 518       if ( currentIndex>=maxhits){
 519         maxhits+=kcmaxhits;
 520         xxx.ResizeTo(4*maxhits);
 521       }
 522       xxx(currentIndex*4)=x[0];
 523       xxx(currentIndex*4+1)=x[1];
 524       xxx(currentIndex*4+2)=x[2];
 525       xxx(currentIndex*4+3)=tpcHit->fQ;
 526       currentIndex++;
 527     }
 528     else
 529       if (currentIndex>2){
 530         Float_t sumx=0;
 531         Float_t sumx2=0;
 532         Float_t sumx3=0;
 533         Float_t sumx4=0;
 534         Float_t sumy=0;
 535         Float_t sumxy=0;
 536         Float_t sumx2y=0;
 537         Float_t sumz=0;
 538         Float_t sumxz=0;
 539         Float_t sumx2z=0;
 540         Float_t sumq=0;
 541         for (Int_t index=0;index<currentIndex;index++){
 542           Float_t x,x2,x3,x4;
 543           x=x2=x3=x4=xxx(index*4);
 544           x2*=x;
 545           x3*=x2;
 546           x4*=x3;
 547           sumx+=x;
 548           sumx2+=x2;
 549           sumx3+=x3;
 550           sumx4+=x4;
 551           sumy+=xxx(index*4+1);
 552           sumxy+=xxx(index*4+1)*x;
 553           sumx2y+=xxx(index*4+1)*x2;
 554           sumz+=xxx(index*4+2);
 555           sumxz+=xxx(index*4+2)*x;
 556           sumx2z+=xxx(index*4+2)*x2;
 557           sumq+=xxx(index*4+3);
 558         }
 559         Float_t centralPad = (param->GetNPads(sector,lastrow)-1)/2;
 560         Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
 561           sumx2*(sumx*sumx3-sumx2*sumx2);
 562
 563         Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
 564           sumx2*(sumxy*sumx3-sumx2y*sumx2);
 565         Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
 566           sumx2*(sumxz*sumx3-sumx2z*sumx2);
 567
 568         Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
 569           sumx2*(sumx*sumx2y-sumx2*sumxy);
 570         Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
 571           sumx2*(sumx*sumx2z-sumx2*sumxz);
 572
 573         Float_t y=detay/det+centralPad;
 574         Float_t z=detaz/det;
 575         Float_t by=detby/det; //y angle
 576         Float_t bz=detbz/det; //z angle
 577         sumy/=Float_t(currentIndex);
 578         sumz/=Float_t(currentIndex);
 579
 580         AliComplexCluster cl;
 581         cl.fX=z;
 582         cl.fY=y;
 583         cl.fQ=sumq;
 584         cl.fSigmaX2=bz;
 585         cl.fSigmaY2=by;
 586         cl.fTracks[0]=ipart;
 587
 588         AliTPCClustersRow * row = (clustersArray->GetRow(sector,lastrow));
 589         if (row!=0) row->InsertCluster(&cl);
 590         currentIndex=0;
 591         lastrow=currentrow;
 592       } //end of calculating cluster for given row
 593
 594   } // end of loop over hits
 595   xxxx->Delete();
 596
 597 }
 598