TPC/AliTPCclustererMI.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 //          Implementation of the TPC clusterer
  20 //
  21 //   Origin: Marian Ivanov
  22 //-------------------------------------------------------
  23
  24 #include "AliTPCReconstructor.h"
  25 #include "AliTPCclustererMI.h"
  26 #include "AliTPCclusterMI.h"
  27 #include <TObjArray.h>
  28 #include <TFile.h>
  29 #include "AliTPCClustersArray.h"
  30 #include "AliTPCClustersRow.h"
  31 #include "AliTPCRawStream.h"
  32 #include "AliDigits.h"
  33 #include "AliSimDigits.h"
  34 #include "AliTPCParam.h"
  35 #include "AliRawReader.h"
  36 #include "AliTPCRawStream.h"
  37 #include "AliRunLoader.h"
  38 #include "AliLoader.h"
  39 #include "Riostream.h"
  40 #include <TTree.h>
  41
  42 ClassImp(AliTPCclustererMI)
  43
  44
  45
  46 AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par)
  47 {
  48   fInput =0;
  49   fOutput=0;
  50   fParam = par;
  51 }
  52 void AliTPCclustererMI::SetInput(TTree * tree)
  53 {
  54   //
  55   // set input tree with digits
  56   //
  57   fInput = tree;
  58   if  (!fInput->GetBranch("Segment")){
  59     cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
  60     fInput=0;
  61     return;
  62   }
  63 }
  64
  65 void AliTPCclustererMI::SetOutput(TTree * tree)
  66 {
  67   //
  68   //
  69   fOutput= tree;
  70   AliTPCClustersRow clrow;
  71   AliTPCClustersRow *pclrow=&clrow;
  72   clrow.SetClass("AliTPCclusterMI");
  73   clrow.SetArray(1); // to make Clones array
  74   fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
  75 }
  76
  77
  78 Float_t  AliTPCclustererMI::GetSigmaY2(Int_t iz){
  79   // sigma y2 = in digits  - we don't know the angle
  80   Float_t z = iz*fParam->GetZWidth();
  81   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
  82     (fPadWidth*fPadWidth);
  83   Float_t sres = 0.25;
  84   Float_t res = sd2+sres;
  85   return res;
  86 }
  87
  88
  89 Float_t  AliTPCclustererMI::GetSigmaZ2(Int_t iz){
  90   //sigma z2 = in digits - angle estimated supposing vertex constraint
  91   Float_t z = iz*fZWidth;
  92   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
  93   Float_t angular = fPadLength*(fParam->GetZLength()-z)/(fRx*fZWidth);
  94   angular*=angular;
  95   angular/=12.;
  96   Float_t sres = fParam->GetZSigma()/fZWidth;
  97   sres *=sres;
  98   Float_t res = angular +sd2+sres;
  99   return res;
 100 }
 101
 102 void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Int_t *bins, UInt_t /*m*/,
 103 AliTPCclusterMI &c)
 104 {
 105   Int_t i0=k/max;  //central pad
 106   Int_t j0=k%max;  //central time bin
 107
 108   // set pointers to data
 109   //Int_t dummy[5] ={0,0,0,0,0};
 110   Int_t * matrix[5]; //5x5 matrix with digits  - indexing i = 0 ..4  j = -2..2
 111   Int_t * resmatrix[5];
 112   for (Int_t di=-2;di<=2;di++){
 113     matrix[di+2]  =  &bins[k+di*max];
 114     resmatrix[di+2]  =  &fResBins[k+di*max];
 115   }
 116   //build matrix with virtual charge
 117   Float_t sigmay2= GetSigmaY2(j0);
 118   Float_t sigmaz2= GetSigmaZ2(j0);
 119
 120   Float_t vmatrix[5][5];
 121   vmatrix[2][2] = matrix[2][0];
 122   c.SetType(0);
 123   c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
 124   for (Int_t di =-1;di <=1;di++)
 125     for (Int_t dj =-1;dj <=1;dj++){
 126       Float_t amp = matrix[di+2][dj];
 127       if ( (amp<2) && (fLoop<2)){
 128         // if under threshold  - calculate virtual charge
 129         Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
 130         amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
 131         if (amp>2) amp = 2;
 132         vmatrix[2+di][2+dj]=amp;
 133         vmatrix[2+2*di][2+2*dj]=0;
 134         if ( (di*dj)!=0){
 135           //DIAGONAL ELEMENTS
 136           vmatrix[2+2*di][2+dj] =0;
 137           vmatrix[2+di][2+2*dj] =0;
 138         }
 139         continue;
 140       }
 141       if (amp<4){
 142         //if small  amplitude - below  2 x threshold  - don't consider other one
 143         vmatrix[2+di][2+dj]=amp;
 144         vmatrix[2+2*di][2+2*dj]=0;  // don't take to the account next bin
 145         if ( (di*dj)!=0){
 146           //DIAGONAL ELEMENTS
 147           vmatrix[2+2*di][2+dj] =0;
 148           vmatrix[2+di][2+2*dj] =0;
 149         }
 150         continue;
 151       }
 152       //if bigger then take everything
 153       vmatrix[2+di][2+dj]=amp;
 154       vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;
 155       if ( (di*dj)!=0){
 156           //DIAGONAL ELEMENTS
 157           vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
 158           vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
 159         }
 160     }
 161
 162
 163
 164   Float_t sumw=0;
 165   Float_t sumiw=0;
 166   Float_t sumi2w=0;
 167   Float_t sumjw=0;
 168   Float_t sumj2w=0;
 169   //
 170   for (Int_t i=-2;i<=2;i++)
 171     for (Int_t j=-2;j<=2;j++){
 172       Float_t amp = vmatrix[i+2][j+2];
 173
 174       sumw    += amp;
 175       sumiw   += i*amp;
 176       sumi2w  += i*i*amp;
 177       sumjw   += j*amp;
 178       sumj2w  += j*j*amp;
 179     }
 180   //
 181   Float_t meani = sumiw/sumw;
 182   Float_t mi2   = sumi2w/sumw-meani*meani;
 183   Float_t meanj = sumjw/sumw;
 184   Float_t mj2   = sumj2w/sumw-meanj*meanj;
 185   //
 186   Float_t ry = mi2/sigmay2;
 187   Float_t rz = mj2/sigmaz2;
 188
 189   //
 190   if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
 191   if ( (ry <1.2) && (rz<1.2) ) {
 192     //if cluster looks like expected
 193     //+1.2 deviation from expected sigma accepted
 194     //    c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
 195
 196     meani +=i0;
 197     meanj +=j0;
 198     //set cluster parameters
 199     c.SetQ(sumw);
 200     c.SetY(meani*fPadWidth);
 201     c.SetZ(meanj*fZWidth);
 202     c.SetSigmaY2(mi2);
 203     c.SetSigmaZ2(mj2);
 204     AddCluster(c);
 205     //remove cluster data from data
 206     for (Int_t di=-2;di<=2;di++)
 207       for (Int_t dj=-2;dj<=2;dj++){
 208         resmatrix[di+2][dj] -= Int_t(vmatrix[di+2][dj+2]);
 209         if (resmatrix[di+2][dj]<0) resmatrix[di+2][dj]=0;
 210       }
 211     resmatrix[2][0] =0;
 212     return;
 213   }
 214   //
 215   //unfolding when neccessary
 216   //
 217
 218   Int_t * matrix2[7]; //7x7 matrix with digits  - indexing i = 0 ..6  j = -3..3
 219   Int_t dummy[7]={0,0,0,0,0,0};
 220   for (Int_t di=-3;di<=3;di++){
 221     matrix2[di+3] =  &bins[k+di*max];
 222     if ((k+di*max)<3)  matrix2[di+3] = &dummy[3];
 223     if ((k+di*max)>fMaxBin-3)  matrix2[di+3] = &dummy[3];
 224   }
 225   Float_t vmatrix2[5][5];
 226   Float_t sumu;
 227   Float_t overlap;
 228   UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
 229   //
 230   //  c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
 231   meani +=i0;
 232   meanj +=j0;
 233   //set cluster parameters
 234   c.SetQ(sumu);
 235   c.SetY(meani*fPadWidth);
 236   c.SetZ(meanj*fZWidth);
 237   c.SetSigmaY2(mi2);
 238   c.SetSigmaZ2(mj2);
 239   c.SetType(Char_t(overlap)+1);
 240   AddCluster(c);
 241
 242   //unfolding 2
 243   meani-=i0;
 244   meanj-=j0;
 245   if (gDebug>4)
 246     printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
 247 }
 248
 249
 250
 251 void AliTPCclustererMI::UnfoldCluster(Int_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
 252                                       Float_t & sumu, Float_t & overlap )
 253 {
 254   //
 255   //unfold cluster from input matrix
 256   //data corresponding to cluster writen in recmatrix
 257   //output meani and meanj
 258
 259   //take separatelly y and z
 260
 261   Float_t sum3i[7] = {0,0,0,0,0,0,0};
 262   Float_t sum3j[7] = {0,0,0,0,0,0,0};
 263
 264   for (Int_t k =0;k<7;k++)
 265     for (Int_t l = -1; l<=1;l++){
 266       sum3i[k]+=matrix2[k][l];
 267       sum3j[k]+=matrix2[l+3][k-3];
 268     }
 269   Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
 270   //
 271   //unfold  y
 272   Float_t sum3wi    = 0;  //charge minus overlap
 273   Float_t sum3wio   = 0;  //full charge
 274   Float_t sum3iw    = 0;  //sum for mean value
 275   for (Int_t dk=-1;dk<=1;dk++){
 276     sum3wio+=sum3i[dk+3];
 277     if (dk==0){
 278       sum3wi+=sum3i[dk+3];
 279     }
 280     else{
 281       Float_t ratio =1;
 282       if (  ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
 283            sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
 284         Float_t xm2 = sum3i[-dk+3];
 285         Float_t xm1 = sum3i[+3];
 286         Float_t x1  = sum3i[2*dk+3];
 287         Float_t x2  = sum3i[3*dk+3];
 288         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
 289         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
 290         ratio = w11/(w11+w12);
 291         for (Int_t dl=-1;dl<=1;dl++)
 292           mratio[dk+1][dl+1] *= ratio;
 293       }
 294       Float_t amp = sum3i[dk+3]*ratio;
 295       sum3wi+=amp;
 296       sum3iw+= dk*amp;
 297     }
 298   }
 299   meani = sum3iw/sum3wi;
 300   Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
 301
 302
 303
 304   //unfold  z
 305   Float_t sum3wj    = 0;  //charge minus overlap
 306   Float_t sum3wjo   = 0;  //full charge
 307   Float_t sum3jw    = 0;  //sum for mean value
 308   for (Int_t dk=-1;dk<=1;dk++){
 309     sum3wjo+=sum3j[dk+3];
 310     if (dk==0){
 311       sum3wj+=sum3j[dk+3];
 312     }
 313     else{
 314       Float_t ratio =1;
 315       if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
 316            (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
 317         Float_t xm2 = sum3j[-dk+3];
 318         Float_t xm1 = sum3j[+3];
 319         Float_t x1  = sum3j[2*dk+3];
 320         Float_t x2  = sum3j[3*dk+3];
 321         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
 322         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
 323         ratio = w11/(w11+w12);
 324         for (Int_t dl=-1;dl<=1;dl++)
 325           mratio[dl+1][dk+1] *= ratio;
 326       }
 327       Float_t amp = sum3j[dk+3]*ratio;
 328       sum3wj+=amp;
 329       sum3jw+= dk*amp;
 330     }
 331   }
 332   meanj = sum3jw/sum3wj;
 333   Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;
 334   overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);
 335   sumu = (sum3wj+sum3wi)/2.;
 336
 337   if (overlap ==3) {
 338     //if not overlap detected remove everything
 339     for (Int_t di =-2; di<=2;di++)
 340       for (Int_t dj =-2; dj<=2;dj++){
 341         recmatrix[di+2][dj+2] = matrix2[3+di][dj];
 342       }
 343   }
 344   else{
 345     for (Int_t di =-1; di<=1;di++)
 346       for (Int_t dj =-1; dj<=1;dj++){
 347         Float_t ratio =1;
 348         if (mratio[di+1][dj+1]==1){
 349           recmatrix[di+2][dj+2]     = matrix2[3+di][dj];
 350           if (TMath::Abs(di)+TMath::Abs(dj)>1){
 351             recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
 352             recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
 353           }
 354           recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
 355         }
 356         else
 357           {
 358             //if we have overlap in direction
 359             recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];
 360             if (TMath::Abs(di)+TMath::Abs(dj)>1){
 361               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
 362               recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
 363               //
 364               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
 365               recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
 366             }
 367             else{
 368               ratio =  recmatrix[di+2][dj+2]/matrix2[3][0];
 369               recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
 370             }
 371           }
 372       }
 373   }
 374   if (gDebug>4)
 375     printf("%f\n", recmatrix[2][2]);
 376
 377 }
 378
 379 Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
 380 {
 381   //
 382   // estimate max
 383   Float_t sumteor= 0;
 384   Float_t sumamp = 0;
 385
 386   for (Int_t di = -1;di<=1;di++)
 387     for (Int_t dj = -1;dj<=1;dj++){
 388       if (vmatrix[2+di][2+dj]>2){
 389         Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
 390         sumteor += teor*vmatrix[2+di][2+dj];
 391         sumamp  += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
 392       }
 393     }
 394   Float_t max = sumamp/sumteor;
 395   return max;
 396 }
 397
 398 void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){
 399   //
 400   // transform cluster to the global coordinata
 401   // add the cluster to the array
 402   //
 403   Float_t meani = c.GetY()/fPadWidth;
 404   Float_t meanj = c.GetZ()/fZWidth;
 405
 406   Int_t ki = TMath::Nint(meani-3);
 407   if (ki<0) ki=0;
 408   if (ki>=fMaxPad) ki = fMaxPad-1;
 409   Int_t kj = TMath::Nint(meanj-3);
 410   if (kj<0) kj=0;
 411   if (kj>=fMaxTime-3) kj=fMaxTime-4;
 412   // ki and kj shifted to "real" coordinata
 413   if (fRowDig) {
 414     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
 415     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
 416     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
 417   }
 418
 419
 420   Float_t s2 = c.GetSigmaY2();
 421   Float_t w=fParam->GetPadPitchWidth(fSector);
 422
 423   c.SetSigmaY2(s2*w*w);
 424   s2 = c.GetSigmaZ2();
 425   w=fZWidth;
 426   c.SetSigmaZ2(s2*w*w);
 427   c.SetY((meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
 428   c.SetZ(fZWidth*(meanj-3));
 429   c.SetZ(c.GetZ()   - 3.*fParam->GetZSigma()); // PASA delay
 430   c.SetZ(fSign*(fParam->GetZLength() - c.GetZ()));
 431
 432   if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
 433     //c.SetSigmaY2(c.GetSigmaY2()*25.);
 434     //c.SetSigmaZ2(c.GetSigmaZ2()*4.);
 435     c.SetType(-(c.GetType()+3));  //edge clusters
 436   }
 437   if (fLoop==2) c.SetType(100);
 438
 439   TClonesArray * arr = fRowCl->GetArray();
 440   // AliTPCclusterMI * cl =
 441   new ((*arr)[fNcluster]) AliTPCclusterMI(c);
 442
 443   fNcluster++;
 444 }
 445
 446
 447 //_____________________________________________________________________________
 448 void AliTPCclustererMI::Digits2Clusters()
 449 {
 450   //-----------------------------------------------------------------
 451   // This is a simple cluster finder.
 452   //-----------------------------------------------------------------
 453
 454   if (!fInput) {
 455     Error("Digits2Clusters", "input tree not initialised");
 456     return;
 457   }
 458
 459   if (!fOutput) {
 460     Error("Digits2Clusters", "output tree not initialised");
 461     return;
 462   }
 463
 464   AliSimDigits digarr, *dummy=&digarr;
 465   fRowDig = dummy;
 466   fInput->GetBranch("Segment")->SetAddress(&dummy);
 467   Stat_t nentries = fInput->GetEntries();
 468
 469   fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3   after
 470
 471   Int_t nclusters  = 0;
 472
 473   for (Int_t n=0; n<nentries; n++) {
 474     fInput->GetEvent(n);
 475     Int_t row;
 476     if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,row)) {
 477       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
 478       continue;
 479     }
 480
 481     AliTPCClustersRow *clrow= new AliTPCClustersRow();
 482     fRowCl = clrow;
 483     clrow->SetClass("AliTPCclusterMI");
 484     clrow->SetArray(1);
 485
 486     clrow->SetID(digarr.GetID());
 487     fOutput->GetBranch("Segment")->SetAddress(&clrow);
 488     fRx=fParam->GetPadRowRadii(fSector,row);
 489
 490
 491     const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
 492     fZWidth = fParam->GetZWidth();
 493     if (fSector < kNIS) {
 494       fMaxPad = fParam->GetNPadsLow(row);
 495       fSign = (fSector < kNIS/2) ? 1 : -1;
 496       fPadLength = fParam->GetPadPitchLength(fSector,row);
 497       fPadWidth = fParam->GetPadPitchWidth();
 498     } else {
 499       fMaxPad = fParam->GetNPadsUp(row);
 500       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 501       fPadLength = fParam->GetPadPitchLength(fSector,row);
 502       fPadWidth  = fParam->GetPadPitchWidth();
 503     }
 504
 505
 506     fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 507     fBins    =new Int_t[fMaxBin];
 508     fResBins =new Int_t[fMaxBin];  //fBins with residuals after 1 finder loop
 509     memset(fBins,0,sizeof(Int_t)*fMaxBin);
 510
 511     if (digarr.First()) //MI change
 512       do {
 513         Short_t dig=digarr.CurrentDigit();
 514         if (dig<=fParam->GetZeroSup()) continue;
 515         Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
 516         fBins[i*fMaxTime+j]=dig;
 517       } while (digarr.Next());
 518     digarr.ExpandTrackBuffer();
 519
 520     FindClusters();
 521
 522     fOutput->Fill();
 523     delete clrow;
 524     nclusters+=fNcluster;
 525     delete[] fBins;
 526     delete[] fResBins;
 527   }
 528
 529   Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
 530 }
 531
 532 void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
 533 {
 534 //-----------------------------------------------------------------
 535 // This is a cluster finder for raw data.
 536 //-----------------------------------------------------------------
 537
 538   if (!fOutput) {
 539     Error("Digits2Clusters", "output tree not initialised");
 540     return;
 541   }
 542
 543   rawReader->Reset();
 544   AliTPCRawStream input(rawReader);
 545
 546   fRowDig = NULL;
 547
 548   Int_t nclusters  = 0;
 549
 550   fMaxTime = fParam->GetMaxTBin() + 6; // add 3 virtual time bins before and 3 after
 551   const Int_t kNIS = fParam->GetNInnerSector();
 552   const Int_t kNOS = fParam->GetNOuterSector();
 553   const Int_t kNS = kNIS + kNOS;
 554   fZWidth = fParam->GetZWidth();
 555   Int_t zeroSup = fParam->GetZeroSup();
 556
 557   fBins = NULL;
 558   Int_t** splitRows = new Int_t* [kNS*2];
 559   Int_t** splitRowsRes = new Int_t* [kNS*2];
 560   for (Int_t iSector = 0; iSector < kNS*2; iSector++)
 561     splitRows[iSector] = NULL;
 562   Int_t iSplitRow = -1;
 563
 564   Bool_t next = kTRUE;
 565   while (next) {
 566     next = input.Next();
 567
 568     // when the sector or row number has changed ...
 569     if (input.IsNewRow() || !next) {
 570
 571       // ... find clusters in the previous pad row, and ...
 572       if (fBins) {
 573         if ((iSplitRow < 0) || splitRows[fSector + kNS*iSplitRow]) {
 574           fRowCl = new AliTPCClustersRow;
 575           fRowCl->SetClass("AliTPCclusterMI");
 576           fRowCl->SetArray(1);
 577           fRowCl->SetID(fParam->GetIndex(fSector, input.GetPrevRow()));
 578           fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
 579
 580           FindClusters();
 581
 582           fOutput->Fill();
 583           delete fRowCl;
 584           nclusters += fNcluster;
 585           delete[] fBins;
 586           delete[] fResBins;
 587           if (iSplitRow >= 0) splitRows[fSector + kNS*iSplitRow] = NULL;
 588
 589         } else if (iSplitRow >= 0) {
 590           splitRows[fSector + kNS*iSplitRow] = fBins;
 591           splitRowsRes[fSector + kNS*iSplitRow] = fResBins;
 592         }
 593       }
 594
 595       if (!next) break;
 596
 597       // ... prepare for the next pad row
 598       fSector = input.GetSector();
 599       Int_t iRow = input.GetRow();
 600       fRx = fParam->GetPadRowRadii(fSector, iRow);
 601
 602       iSplitRow = -1;
 603       if (fSector < kNIS) {
 604         fMaxPad = fParam->GetNPadsLow(iRow);
 605         fSign = (fSector < kNIS/2) ? 1 : -1;
 606         if (iRow == 30) iSplitRow = 0;
 607       } else {
 608         fMaxPad = fParam->GetNPadsUp(iRow);
 609         fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 610         if (iRow == 27) iSplitRow = 0;
 611         else if (iRow == 76) iSplitRow = 1;
 612       }
 613       fPadLength = fParam->GetPadPitchLength(fSector, iRow);
 614       fPadWidth  = fParam->GetPadPitchWidth();
 615
 616       fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 617       if ((iSplitRow < 0) || !splitRows[fSector + kNS*iSplitRow]) {
 618         fBins    = new Int_t[fMaxBin];
 619         fResBins = new Int_t[fMaxBin];  //fBins with residuals after 1 finder loop
 620         memset(fBins, 0, sizeof(Int_t)*fMaxBin);
 621       } else {
 622         fBins    = splitRows[fSector + kNS*iSplitRow];
 623         fResBins = splitRowsRes[fSector + kNS*iSplitRow];
 624       }
 625     }
 626
 627     // fill fBins with digits data
 628     if (input.GetSignal() <= zeroSup) continue;
 629     Int_t i = input.GetPad() + 3;
 630     Int_t j = input.GetTime() + 3;
 631     fBins[i*fMaxTime+j] = input.GetSignal();
 632   }
 633
 634   // find clusters in split rows that were skipped until now.
 635   // this can happen if the rows were not splitted
 636   for (fSector = 0; fSector < kNS; fSector++)
 637     for (Int_t iSplit = 0; iSplit < 2; iSplit++)
 638       if (splitRows[fSector + kNS*iSplit]) {
 639
 640         Int_t iRow = -1;
 641         if (fSector < kNIS) {
 642           iRow = 30;
 643           fMaxPad = fParam->GetNPadsLow(iRow);
 644           fSign = (fSector < kNIS/2) ? 1 : -1;
 645         } else {
 646           if (iSplit == 0) iRow = 27; else iRow = 76;
 647           fMaxPad = fParam->GetNPadsUp(iRow);
 648           fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 649         }
 650         fRx = fParam->GetPadRowRadii(fSector, iRow);
 651         fPadLength = fParam->GetPadPitchLength(fSector, iRow);
 652         fPadWidth  = fParam->GetPadPitchWidth();
 653
 654         fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 655         fBins    = splitRows[fSector + kNS*iSplit];
 656         fResBins = splitRowsRes[fSector + kNS*iSplit];
 657
 658         fRowCl = new AliTPCClustersRow;
 659         fRowCl->SetClass("AliTPCclusterMI");
 660         fRowCl->SetArray(1);
 661         fRowCl->SetID(fParam->GetIndex(fSector, iRow));
 662         fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
 663
 664         FindClusters();
 665
 666         fOutput->Fill();
 667         delete fRowCl;
 668         nclusters += fNcluster;
 669         delete[] fBins;
 670         delete[] fResBins;
 671       }
 672
 673   delete[] splitRows;
 674   delete[] splitRowsRes;
 675   Info("Digits2Clusters", "Number of found clusters : %d\n", nclusters);
 676 }
 677
 678 void AliTPCclustererMI::FindClusters()
 679 {
 680   //add virtual charge at the edge
 681   for (Int_t i=0; i<fMaxTime; i++){
 682     Float_t amp1 = fBins[i+3*fMaxTime];
 683     Float_t amp0 =0;
 684     if (amp1>0){
 685       Float_t amp2 = fBins[i+4*fMaxTime];
 686       if (amp2==0) amp2=0.5;
 687       Float_t sigma2 = GetSigmaY2(i);
 688       amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);
 689       if (gDebug>4) printf("\n%f\n",amp0);
 690     }
 691     fBins[i+2*fMaxTime] = Int_t(amp0);
 692     amp0 = 0;
 693     amp1 = fBins[(fMaxPad+2)*fMaxTime+i];
 694     if (amp1>0){
 695       Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime];
 696       if (amp2==0) amp2=0.5;
 697       Float_t sigma2 = GetSigmaY2(i);
 698       amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);
 699       if (gDebug>4) printf("\n%f\n",amp0);
 700     }
 701     fBins[(fMaxPad+3)*fMaxTime+i] = Int_t(amp0);
 702   }
 703
 704 //  memcpy(fResBins,fBins, fMaxBin*2);
 705   memcpy(fResBins,fBins, fMaxBin);
 706   //
 707   fNcluster=0;
 708   //first loop - for "gold cluster"
 709   fLoop=1;
 710   Int_t *b=&fBins[-1]+2*fMaxTime;
 711   Int_t crtime = Int_t((fParam->GetZLength()-AliTPCReconstructor::GetCtgRange()*fRx)/fZWidth-5);
 712
 713   for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
 714     b++;
 715     if (*b<8) continue;   //threshold form maxima
 716     if (i%fMaxTime<crtime) {
 717       Int_t delta = -(i%fMaxTime)+crtime;
 718       b+=delta;
 719       i+=delta;
 720       continue;
 721     }
 722
 723     if (!IsMaximum(*b,fMaxTime,b)) continue;
 724     AliTPCclusterMI c;
 725     Int_t dummy=0;
 726     MakeCluster(i, fMaxTime, fBins, dummy,c);
 727     //}
 728   }
 729   //memcpy(fBins,fResBins, fMaxBin*2);
 730   //second  loop - for rest cluster
 731   /*
 732   fLoop=2;
 733   b=&fResBins[-1]+2*fMaxTime;
 734   for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
 735     b++;
 736     if (*b<25) continue;   // bigger threshold for maxima
 737     if (!IsMaximum(*b,fMaxTime,b)) continue;
 738     AliTPCclusterMI c;
 739     Int_t dummy;
 740     MakeCluster(i, fMaxTime, fResBins, dummy,c);
 741     //}
 742   }
 743   */
 744 }