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 //  1. The Input data for reconstruction - Options
  22 //      1.a Simulated data  - TTree - invoked Digits2Clusters()
  23 //      1.b Raw data        - Digits2Clusters(AliRawReader* rawReader);
  24 //      1.c HLT clusters    - Digits2Clusters and Digits2Clusters(AliRawReader* rawReader)
  25 //                            invoke ReadHLTClusters()
  26 //
  27 //      fUseHLTClusters     - switches between different inputs
  28 //                            1 -> only TPC raw/sim data
  29 //                            2 -> if present TPC raw/sim data, otherwise HLT clusters
  30 //                            3 -> only HLT clusters
  31 //                            4 -> if present HLT clusters, otherwise TPC raw/sim data
  32 //
  33 //  2. The Output data
  34 //      2.a TTree with clusters - if  SetOutput(TTree * tree) invoked
  35 //      2.b TObjArray           - Faster option for HLT
  36 //      2.c TClonesArray        - Faster option for HLT (smaller memory consumption), activate with fBClonesArray flag
  37 //
  38 //  3. Reconstruction setup
  39 //     see AliTPCRecoParam for list of parameters
  40 //     The reconstruction parameterization taken from the
  41 //     AliTPCReconstructor::GetRecoParam()
  42 //     Possible to setup it in reconstruction macro  AliTPCReconstructor::SetRecoParam(...)
  43 //
  44 //
  45 //
  46 //   Origin: Marian Ivanov
  47 //-------------------------------------------------------
  48
  49 #include "Riostream.h"
  50 #include <TF1.h>
  51 #include <TFile.h>
  52 #include <TGraph.h>
  53 #include <TH1F.h>
  54 #include <TObjArray.h>
  55 #include <TClonesArray.h>
  56 #include <TRandom.h>
  57 #include <TTree.h>
  58 #include <TTreeStream.h>
  59 #include "TSystem.h"
  60 #include "TClass.h"
  61
  62 #include "AliDigits.h"
  63 #include "AliLoader.h"
  64 #include "AliLog.h"
  65 #include "AliMathBase.h"
  66 #include "AliRawEventHeaderBase.h"
  67 #include "AliRawReader.h"
  68 #include "AliRunLoader.h"
  69 #include "AliSimDigits.h"
  70 #include "AliTPCCalPad.h"
  71 #include "AliTPCCalROC.h"
  72 #include "AliTPCClustersArray.h"
  73 #include "AliTPCClustersRow.h"
  74 #include "AliTPCParam.h"
  75 #include "AliTPCRawStream.h"
  76 #include "AliTPCRawStreamV3.h"
  77 #include "AliTPCRecoParam.h"
  78 #include "AliTPCReconstructor.h"
  79 #include "AliTPCcalibDB.h"
  80 #include "AliTPCclusterInfo.h"
  81 #include "AliTPCclusterMI.h"
  82 #include "AliTPCTransform.h"
  83 #include "AliTPCclustererMI.h"
  84
  85 using std::cerr;
  86 using std::endl;
  87 ClassImp(AliTPCclustererMI)
  88
  89
  90
  91 AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam):
  92   fBins(0),
  93   fSigBins(0),
  94   fNSigBins(0),
  95   fLoop(0),
  96   fMaxBin(0),
  97   fMaxTime(1006), // 1000>940 so use 1000, add 3 virtual time bins before and 3 after
  98   fMaxPad(0),
  99   fSector(-1),
 100   fRow(-1),
 101   fSign(0),
 102   fRx(0),
 103   fPadWidth(0),
 104   fPadLength(0),
 105   fZWidth(0),
 106   fPedSubtraction(kFALSE),
 107   fEventHeader(0),
 108   fTimeStamp(0),
 109   fEventType(0),
 110   fInput(0),
 111   fOutput(0),
 112   fOutputArray(0),
 113   fOutputClonesArray(0),
 114   fRowCl(0),
 115   fRowDig(0),
 116   fParam(0),
 117   fNcluster(0),
 118   fNclusters(0),
 119   fDebugStreamer(0),
 120   fRecoParam(0),
 121   fBDumpSignal(kFALSE),
 122   fBClonesArray(kFALSE),
 123   fUseHLTClusters(4),
 124   fAllBins(NULL),
 125   fAllSigBins(NULL),
 126   fAllNSigBins(NULL),
 127   fHLTClusterAccess(NULL)
 128 {
 129   //
 130   // COSNTRUCTOR
 131   // param     - tpc parameters for given file
 132   // recoparam - reconstruction parameters
 133   //
 134   fInput =0;
 135   fParam = par;
 136   if (recoParam) {
 137     fRecoParam = recoParam;
 138   }else{
 139     //set default parameters if not specified
 140     fRecoParam = AliTPCReconstructor::GetRecoParam();
 141     if (!fRecoParam)  fRecoParam = AliTPCRecoParam::GetLowFluxParam();
 142   }
 143
 144   if(AliTPCReconstructor::StreamLevel()>0) {
 145     fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
 146   }
 147
 148   //  Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
 149   fRowCl= new AliTPCClustersRow("AliTPCclusterMI");
 150
 151   // Non-persistent arrays
 152   //
 153   //alocate memory for sector - maximal case
 154   //
 155   AliTPCROC * roc = AliTPCROC::Instance();
 156   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
 157   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
 158
 159   fAllBins = new Float_t*[nRowsMax];
 160   fAllSigBins = new Int_t*[nRowsMax];
 161   fAllNSigBins = new Int_t[nRowsMax];
 162   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
 163     //
 164     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
 165     fAllBins[iRow] = new Float_t[maxBin];
 166     memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
 167     fAllSigBins[iRow] = new Int_t[maxBin];
 168     fAllNSigBins[iRow]=0;
 169   }
 170 }
 171
 172 //______________________________________________________________
 173 AliTPCclustererMI::~AliTPCclustererMI(){
 174   //
 175   //
 176   //
 177   if (fDebugStreamer) delete fDebugStreamer;
 178   if (fOutputArray){
 179     //fOutputArray->Delete();
 180     delete fOutputArray;
 181   }
 182   if (fOutputClonesArray){
 183     fOutputClonesArray->Delete();
 184     delete fOutputClonesArray;
 185   }
 186
 187   if (fRowCl) {
 188     fRowCl->GetArray()->Delete();
 189     delete fRowCl;
 190   }
 191
 192   AliTPCROC * roc = AliTPCROC::Instance();
 193   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
 194   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
 195     delete [] fAllBins[iRow];
 196     delete [] fAllSigBins[iRow];
 197   }
 198   delete [] fAllBins;
 199   delete [] fAllSigBins;
 200   delete [] fAllNSigBins;
 201   if (fHLTClusterAccess) delete fHLTClusterAccess;
 202 }
 203
 204 void AliTPCclustererMI::SetInput(TTree * tree)
 205 {
 206   //
 207   // set input tree with digits
 208   //
 209   fInput = tree;
 210   if  (!fInput->GetBranch("Segment")){
 211     cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
 212     fInput=0;
 213     return;
 214   }
 215 }
 216
 217 void AliTPCclustererMI::SetOutput(TTree * tree)
 218 {
 219   //
 220   // Set the output tree
 221   // If not set the ObjArray used - Option for HLT
 222   //
 223   if (!tree) return;
 224   fOutput= tree;
 225   AliTPCClustersRow clrow("AliTPCclusterMI");
 226   AliTPCClustersRow *pclrow=&clrow;
 227   fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
 228 }
 229
 230
 231 void AliTPCclustererMI::FillRow(){
 232   //
 233   // fill the output container -
 234   // 2 Options possible
 235   //          Tree
 236   //          TObjArray
 237   //
 238   if (fOutput) fOutput->Fill();
 239   if (!fOutput && !fBClonesArray){
 240     //
 241     if (!fOutputArray) fOutputArray = new TObjArray(fParam->GetNRowsTotal());
 242     if (fRowCl && fRowCl->GetArray()->GetEntriesFast()>0) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID());
 243   }
 244 }
 245
 246 Float_t  AliTPCclustererMI::GetSigmaY2(Int_t iz){
 247   // sigma y2 = in digits  - we don't know the angle
 248   Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
 249   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
 250     (fPadWidth*fPadWidth);
 251   Float_t sres = 0.25;
 252   Float_t res = sd2+sres;
 253   return res;
 254 }
 255
 256
 257 Float_t  AliTPCclustererMI::GetSigmaZ2(Int_t iz){
 258   //sigma z2 = in digits - angle estimated supposing vertex constraint
 259   Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
 260   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
 261   Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
 262   angular*=angular;
 263   angular/=12.;
 264   Float_t sres = fParam->GetZSigma()/fZWidth;
 265   sres *=sres;
 266   Float_t res = angular +sd2+sres;
 267   return res;
 268 }
 269
 270 void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
 271 AliTPCclusterMI &c)
 272 {
 273   //
 274   //  Make cluster: characterized by position ( mean-  COG) , shape (RMS) a charge, QMax and Q tot
 275   //  Additional correction:
 276   //  a) To correct for charge below threshold, in the +1 neghborhood to the max charge charge
 277   //       is extrapolated using gaussian approximation assuming given cluster width..
 278   //       Additional empirical factor is used to account for the charge fluctuation (kVirtualChargeFactor).
 279   //       Actual value of the  kVirtualChargeFactor should obtained minimimizing residuals between the cluster
 280   //       and track interpolation.
 281   //  b.) For space points with extended shape (in comparison with expected using parameterization) clusters are
 282   //      unfoded
 283   //
 284   //  NOTE. Actual/Empirical  values for correction are hardwired in the code.
 285   //
 286   // Input paramters for function:
 287   //  k    - Make cluster at position k
 288   //  bins - 2 D array of signals mapped to 1 dimensional array -
 289   //  max  - the number of time bins er one dimension
 290   //  c    - reference to cluster to be filled
 291   //
 292   Double_t kVirtualChargeFactor=0.5;
 293   Int_t i0=k/max;  //central pad
 294   Int_t j0=k%max;  //central time bin
 295
 296   // set pointers to data
 297   //Int_t dummy[5] ={0,0,0,0,0};
 298   Float_t * matrix[5]; //5x5 matrix with digits  - indexing i = 0 ..4  j = -2..2
 299   for (Int_t di=-2;di<=2;di++){
 300     matrix[di+2]  =  &bins[k+di*max];
 301   }
 302   //build matrix with virtual charge
 303   Float_t sigmay2= GetSigmaY2(j0);
 304   Float_t sigmaz2= GetSigmaZ2(j0);
 305
 306   Float_t vmatrix[5][5];
 307   vmatrix[2][2] = matrix[2][0];
 308   c.SetType(0);
 309   c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
 310   for (Int_t di =-1;di <=1;di++)
 311     for (Int_t dj =-1;dj <=1;dj++){
 312       Float_t amp = matrix[di+2][dj];
 313       if ( (amp<2) && (fLoop<2)){
 314         // if under threshold  - calculate virtual charge
 315         Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
 316         amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
 317         if (amp>2) amp = 2;
 318         vmatrix[2+di][2+dj]= kVirtualChargeFactor*amp;
 319         vmatrix[2+2*di][2+2*dj]=0;
 320         if ( (di*dj)!=0){
 321           //DIAGONAL ELEMENTS
 322           vmatrix[2+2*di][2+dj] =0;
 323           vmatrix[2+di][2+2*dj] =0;
 324         }
 325         continue;
 326       }
 327       if (amp<4){
 328         //if small  amplitude - below  2 x threshold  - don't consider other one
 329         vmatrix[2+di][2+dj]=amp;
 330         vmatrix[2+2*di][2+2*dj]=0;  // don't take to the account next bin
 331         if ( (di*dj)!=0){
 332           //DIAGONAL ELEMENTS
 333           vmatrix[2+2*di][2+dj] =0;
 334           vmatrix[2+di][2+2*dj] =0;
 335         }
 336         continue;
 337       }
 338       //if bigger then take everything
 339       vmatrix[2+di][2+dj]=amp;
 340       vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;
 341       if ( (di*dj)!=0){
 342           //DIAGONAL ELEMENTS
 343           vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
 344           vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
 345         }
 346     }
 347
 348
 349
 350   Float_t sumw=0;
 351   Float_t sumiw=0;
 352   Float_t sumi2w=0;
 353   Float_t sumjw=0;
 354   Float_t sumj2w=0;
 355   //
 356   for (Int_t i=-2;i<=2;i++)
 357     for (Int_t j=-2;j<=2;j++){
 358       Float_t amp = vmatrix[i+2][j+2];
 359
 360       sumw    += amp;
 361       sumiw   += i*amp;
 362       sumi2w  += i*i*amp;
 363       sumjw   += j*amp;
 364       sumj2w  += j*j*amp;
 365     }
 366   //
 367   Float_t meani = sumiw/sumw;
 368   Float_t mi2   = sumi2w/sumw-meani*meani;
 369   Float_t meanj = sumjw/sumw;
 370   Float_t mj2   = sumj2w/sumw-meanj*meanj;
 371   //
 372   Float_t ry = mi2/sigmay2;
 373   Float_t rz = mj2/sigmaz2;
 374
 375   //
 376   if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
 377   if ( ((ry <1.2) && (rz<1.2)) || (!fRecoParam->GetDoUnfold())) {
 378     //
 379     //if cluster looks like expected or Unfolding not switched on
 380     //standard COG is used
 381     //+1.2 deviation from expected sigma accepted
 382     //    c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
 383
 384     meani +=i0;
 385     meanj +=j0;
 386     //set cluster parameters
 387     c.SetQ(sumw);
 388     c.SetPad(meani-2.5);
 389     c.SetTimeBin(meanj-3);
 390     c.SetSigmaY2(mi2);
 391     c.SetSigmaZ2(mj2);
 392     c.SetType(0);
 393     AddCluster(c,(Float_t*)vmatrix,k);
 394     return;
 395   }
 396   //
 397   //unfolding when neccessary
 398   //
 399
 400   Float_t * matrix2[7]; //7x7 matrix with digits  - indexing i = 0 ..6  j = -3..3
 401   Float_t dummy[7]={0,0,0,0,0,0};
 402   for (Int_t di=-3;di<=3;di++){
 403     matrix2[di+3] =  &bins[k+di*max];
 404     if ((k+di*max)<3)  matrix2[di+3] = &dummy[3];
 405     if ((k+di*max)>fMaxBin-3)  matrix2[di+3] = &dummy[3];
 406   }
 407   Float_t vmatrix2[5][5];
 408   Float_t sumu;
 409   Float_t overlap;
 410   UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
 411   //
 412   //  c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
 413   meani +=i0;
 414   meanj +=j0;
 415   //set cluster parameters
 416   c.SetQ(sumu);
 417   c.SetPad(meani-2.5);
 418   c.SetTimeBin(meanj-3);
 419   c.SetSigmaY2(mi2);
 420   c.SetSigmaZ2(mj2);
 421   c.SetType(Char_t(overlap)+1);
 422   AddCluster(c,(Float_t*)vmatrix,k);
 423
 424   //unfolding 2
 425   meani-=i0;
 426   meanj-=j0;
 427 }
 428
 429
 430
 431 void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
 432                                       Float_t & sumu, Float_t & overlap )
 433 {
 434   //
 435   //unfold cluster from input matrix
 436   //data corresponding to cluster writen in recmatrix
 437   //output meani and meanj
 438
 439   //take separatelly y and z
 440
 441   Float_t sum3i[7] = {0,0,0,0,0,0,0};
 442   Float_t sum3j[7] = {0,0,0,0,0,0,0};
 443
 444   for (Int_t k =0;k<7;k++)
 445     for (Int_t l = -1; l<=1;l++){
 446       sum3i[k]+=matrix2[k][l];
 447       sum3j[k]+=matrix2[l+3][k-3];
 448     }
 449   Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
 450   //
 451   //unfold  y
 452   Float_t sum3wi    = 0;  //charge minus overlap
 453   Float_t sum3wio   = 0;  //full charge
 454   Float_t sum3iw    = 0;  //sum for mean value
 455   for (Int_t dk=-1;dk<=1;dk++){
 456     sum3wio+=sum3i[dk+3];
 457     if (dk==0){
 458       sum3wi+=sum3i[dk+3];
 459     }
 460     else{
 461       Float_t ratio =1;
 462       if (  ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
 463             (sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 )){
 464         Float_t xm2 = sum3i[-dk+3];
 465         Float_t xm1 = sum3i[+3];
 466         Float_t x1  = sum3i[2*dk+3];
 467         Float_t x2  = sum3i[3*dk+3];
 468         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
 469         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
 470         ratio = w11/(w11+w12);
 471         for (Int_t dl=-1;dl<=1;dl++)
 472           mratio[dk+1][dl+1] *= ratio;
 473       }
 474       Float_t amp = sum3i[dk+3]*ratio;
 475       sum3wi+=amp;
 476       sum3iw+= dk*amp;
 477     }
 478   }
 479   meani = sum3iw/sum3wi;
 480   Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
 481
 482
 483
 484   //unfold  z
 485   Float_t sum3wj    = 0;  //charge minus overlap
 486   Float_t sum3wjo   = 0;  //full charge
 487   Float_t sum3jw    = 0;  //sum for mean value
 488   for (Int_t dk=-1;dk<=1;dk++){
 489     sum3wjo+=sum3j[dk+3];
 490     if (dk==0){
 491       sum3wj+=sum3j[dk+3];
 492     }
 493     else{
 494       Float_t ratio =1;
 495       if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
 496            (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
 497         Float_t xm2 = sum3j[-dk+3];
 498         Float_t xm1 = sum3j[+3];
 499         Float_t x1  = sum3j[2*dk+3];
 500         Float_t x2  = sum3j[3*dk+3];
 501         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
 502         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
 503         ratio = w11/(w11+w12);
 504         for (Int_t dl=-1;dl<=1;dl++)
 505           mratio[dl+1][dk+1] *= ratio;
 506       }
 507       Float_t amp = sum3j[dk+3]*ratio;
 508       sum3wj+=amp;
 509       sum3jw+= dk*amp;
 510     }
 511   }
 512   meanj = sum3jw/sum3wj;
 513   Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;
 514   overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);
 515   sumu = (sum3wj+sum3wi)/2.;
 516
 517   if (overlap ==3) {
 518     //if not overlap detected remove everything
 519     for (Int_t di =-2; di<=2;di++)
 520       for (Int_t dj =-2; dj<=2;dj++){
 521         recmatrix[di+2][dj+2] = matrix2[3+di][dj];
 522       }
 523   }
 524   else{
 525     for (Int_t di =-1; di<=1;di++)
 526       for (Int_t dj =-1; dj<=1;dj++){
 527         Float_t ratio =1;
 528         if (mratio[di+1][dj+1]==1){
 529           recmatrix[di+2][dj+2]     = matrix2[3+di][dj];
 530           if (TMath::Abs(di)+TMath::Abs(dj)>1){
 531             recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
 532             recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
 533           }
 534           recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
 535         }
 536         else
 537           {
 538             //if we have overlap in direction
 539             recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];
 540             if (TMath::Abs(di)+TMath::Abs(dj)>1){
 541               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
 542               recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
 543               //
 544               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
 545               recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
 546             }
 547             else{
 548               ratio =  recmatrix[di+2][dj+2]/matrix2[3][0];
 549               recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
 550             }
 551           }
 552       }
 553   }
 554
 555 }
 556
 557 Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
 558 {
 559   //
 560   // estimate max
 561   Float_t sumteor= 0;
 562   Float_t sumamp = 0;
 563
 564   for (Int_t di = -1;di<=1;di++)
 565     for (Int_t dj = -1;dj<=1;dj++){
 566       if (vmatrix[2+di][2+dj]>2){
 567         Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
 568         sumteor += teor*vmatrix[2+di][2+dj];
 569         sumamp  += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
 570       }
 571     }
 572   Float_t max = sumamp/sumteor;
 573   return max;
 574 }
 575
 576 void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * /*matrix*/, Int_t /*pos*/){
 577   //
 578   //
 579   // Transform cluster to the rotated global coordinata
 580   // Assign labels to the cluster
 581   // add the cluster to the array
 582   // for more details - See  AliTPCTranform::Transform(x,i,0,1)
 583   Float_t meani = c.GetPad();
 584   Float_t meanj = c.GetTimeBin();
 585
 586   Int_t ki = TMath::Nint(meani);
 587   if (ki<0) ki=0;
 588   if (ki>=fMaxPad) ki = fMaxPad-1;
 589   Int_t kj = TMath::Nint(meanj);
 590   if (kj<0) kj=0;
 591   if (kj>=fMaxTime-3) kj=fMaxTime-4;
 592   // ki and kj shifted as integers coordinata
 593   if (fRowDig) {
 594     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
 595     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
 596     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
 597   }
 598
 599   c.SetRow(fRow);
 600   c.SetDetector(fSector);
 601   Float_t s2 = c.GetSigmaY2();
 602   Float_t w=fParam->GetPadPitchWidth(fSector);
 603   c.SetSigmaY2(s2*w*w);
 604   s2 = c.GetSigmaZ2();
 605   c.SetSigmaZ2(s2*fZWidth*fZWidth);
 606   //
 607   //
 608   //
 609
 610   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
 611   if (!transform) {
 612     AliFatal("Tranformations not in calibDB");
 613     return;
 614   }
 615   transform->SetCurrentRecoParam((AliTPCRecoParam*)fRecoParam);
 616   Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()};
 617   Int_t i[1]={fSector};
 618   transform->Transform(x,i,0,1);
 619   c.SetX(x[0]);
 620   c.SetY(x[1]);
 621   c.SetZ(x[2]);
 622   //
 623   //
 624   if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
 625     c.SetType(-(c.GetType()+3));  //edge clusters
 626   }
 627   if (fLoop==2) c.SetType(100);
 628
 629   // select output
 630   TClonesArray * arr = 0;
 631   AliTPCclusterMI * cl = 0;
 632
 633   if(fBClonesArray==kFALSE) {
 634      arr = fRowCl->GetArray();
 635      cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
 636   } else {
 637      cl = new ((*fOutputClonesArray)[fNclusters+fNcluster]) AliTPCclusterMI(c);
 638   }
 639
 640   // if (fRecoParam->DumpSignal() &&matrix ) {
 641 //     Int_t nbins=0;
 642 //     Float_t *graph =0;
 643 //     if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
 644 //       nbins = fMaxTime;
 645 //       graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
 646 //     }
 647 //     AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
 648 //     cl->SetInfo(info);
 649 //   }
 650   if (!fRecoParam->DumpSignal()) {
 651     cl->SetInfo(0);
 652   }
 653
 654   if (AliTPCReconstructor::StreamLevel()>1) {
 655     Float_t xyz[3];
 656     cl->GetGlobalXYZ(xyz);
 657      (*fDebugStreamer)<<"Clusters"<<
 658        "Cl.="<<cl<<
 659        "gx="<<xyz[0]<<
 660        "gy="<<xyz[1]<<
 661        "gz="<<xyz[2]<<
 662        "\n";
 663   }
 664
 665   fNcluster++;
 666 }
 667
 668
 669 //_____________________________________________________________________________
 670 void AliTPCclustererMI::Digits2Clusters()
 671 {
 672   //-----------------------------------------------------------------
 673   // This is a simple cluster finder.
 674   //-----------------------------------------------------------------
 675
 676   if (!fInput) {
 677     Error("Digits2Clusters", "input tree not initialised");
 678     return;
 679   }
 680   fRecoParam = AliTPCReconstructor::GetRecoParam();
 681   if (!fRecoParam){
 682     AliFatal("Can not get the reconstruction parameters");
 683   }
 684   if(AliTPCReconstructor::StreamLevel()>5) {
 685     AliInfo("Parameter Dumps");
 686     fParam->Dump();
 687     fRecoParam->Dump();
 688   }
 689   fRowDig = NULL;
 690
 691   //-----------------------------------------------------------------
 692   // Use HLT clusters
 693   //-----------------------------------------------------------------
 694   if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
 695     AliInfo("Using HLT clusters for TPC off-line reconstruction");
 696     fZWidth = fParam->GetZWidth();
 697     Int_t iResult = ReadHLTClusters();
 698
 699     // HLT clusters present
 700     if (iResult >= 0 && fNclusters > 0)
 701       return;
 702
 703     // HLT clusters not present
 704     if (iResult < 0 || fNclusters == 0) {
 705       if (fUseHLTClusters == 3) {
 706         AliError("No HLT clusters present, but requested.");
 707         return;
 708       }
 709       else {
 710         AliInfo("Now trying to read from TPC RAW");
 711       }
 712     }
 713     // Some other problem during cluster reading
 714     else {
 715       AliWarning("Some problem while unpacking of HLT clusters.");
 716       return;
 717     }
 718   } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
 719
 720   //-----------------------------------------------------------------
 721   // Run TPC off-line clusterer
 722   //-----------------------------------------------------------------
 723   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
 724   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 725   AliSimDigits digarr, *dummy=&digarr;
 726   fRowDig = dummy;
 727   fInput->GetBranch("Segment")->SetAddress(&dummy);
 728   Stat_t nentries = fInput->GetEntries();
 729
 730   fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3   after
 731
 732   Int_t nclusters  = 0;
 733
 734   for (Int_t n=0; n<nentries; n++) {
 735     fInput->GetEvent(n);
 736     if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
 737       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
 738       continue;
 739     }
 740     Int_t row = fRow;
 741     AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector);  // pad gains per given sector
 742     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector); // noise per given sector
 743     //
 744
 745     fRowCl->SetID(digarr.GetID());
 746     if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
 747     fRx=fParam->GetPadRowRadii(fSector,row);
 748
 749
 750     const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
 751     fZWidth = fParam->GetZWidth();
 752     if (fSector < kNIS) {
 753       fMaxPad = fParam->GetNPadsLow(row);
 754       fSign =  (fSector < kNIS/2) ? 1 : -1;
 755       fPadLength = fParam->GetPadPitchLength(fSector,row);
 756       fPadWidth = fParam->GetPadPitchWidth();
 757     } else {
 758       fMaxPad = fParam->GetNPadsUp(row);
 759       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 760       fPadLength = fParam->GetPadPitchLength(fSector,row);
 761       fPadWidth  = fParam->GetPadPitchWidth();
 762     }
 763
 764
 765     fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 766     fBins    =new Float_t[fMaxBin];
 767     fSigBins =new Int_t[fMaxBin];
 768     fNSigBins = 0;
 769     memset(fBins,0,sizeof(Float_t)*fMaxBin);
 770
 771     if (digarr.First()) //MI change
 772       do {
 773         Float_t dig=digarr.CurrentDigit();
 774         if (dig<=fParam->GetZeroSup()) continue;
 775         Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
 776         Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
 777         Int_t bin = i*fMaxTime+j;
 778         if (gain>0){
 779           fBins[bin]=dig/gain;
 780         }else{
 781           fBins[bin]=0;
 782         }
 783         fSigBins[fNSigBins++]=bin;
 784       } while (digarr.Next());
 785     digarr.ExpandTrackBuffer();
 786
 787     FindClusters(noiseROC);
 788     FillRow();
 789     fRowCl->GetArray()->Clear("C");
 790     nclusters+=fNcluster;
 791
 792     delete[] fBins;
 793     delete[] fSigBins;
 794   }
 795
 796   Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
 797
 798   if (fUseHLTClusters == 2 && nclusters == 0) {
 799     AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters.");
 800
 801     fZWidth = fParam->GetZWidth();
 802     ReadHLTClusters();
 803   }
 804 }
 805
 806 void AliTPCclustererMI::ProcessSectorData(){
 807   //
 808   // Process the data for the current sector
 809   //
 810
 811   AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
 812   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 813   AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector);  // pedestal per given sector
 814   AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector);  // noise per given sector
 815   //check the presence of the calibration
 816   if (!noiseROC ||!pedestalROC ) {
 817     AliError(Form("Missing calibration per sector\t%d\n",fSector));
 818     return;
 819   }
 820   Int_t  nRows=fParam->GetNRow(fSector);
 821   Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
 822   Int_t zeroSup = fParam->GetZeroSup();
 823   //    if (calcPedestal) {
 824   if (kFALSE ) {
 825     for (Int_t iRow = 0; iRow < nRows; iRow++) {
 826       Int_t maxPad = fParam->GetNPads(fSector, iRow);
 827
 828       for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
 829     //
 830     // Temporary fix for data production - !!!! MARIAN
 831     // The noise calibration should take mean and RMS - currently the Gaussian fit used
 832     // In case of double peak  - the pad should be rejected
 833     //
 834     // Line mean - if more than given digits over threshold - make a noise calculation
 835     // and pedestal substration
 836         if (!calcPedestal && fAllBins[iRow][iPad*fMaxTime+0]<50) continue;
 837     //
 838         if (fAllBins[iRow][iPad*fMaxTime+0] <1 ) continue;  // no data
 839         Float_t *p = &fAllBins[iRow][iPad*fMaxTime+3];
 840     //Float_t pedestal = TMath::Median(fMaxTime, p);
 841         Int_t id[3] = {fSector, iRow, iPad-3};
 842     // calib values
 843         Double_t rmsCalib=  noiseROC->GetValue(iRow,iPad-3);
 844         Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
 845         Double_t rmsEvent       = rmsCalib;
 846         Double_t pedestalEvent  = pedestalCalib;
 847         ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
 848         if (rmsEvent<rmsCalib) rmsEvent = rmsCalib;   // take worst scenario
 849         if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
 850
 851     //
 852         for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
 853           Int_t bin = iPad*fMaxTime+iTimeBin;
 854           fAllBins[iRow][bin] -= pedestalEvent;
 855           if (iTimeBin < fRecoParam->GetFirstBin())
 856             fAllBins[iRow][bin] = 0;
 857           if (iTimeBin > fRecoParam->GetLastBin())
 858             fAllBins[iRow][bin] = 0;
 859           if (fAllBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
 860             fAllBins[iRow][bin] = 0;
 861           if (fAllBins[iRow][bin] < 3.0*rmsEvent)   // 3 sigma cut on RMS
 862             fAllBins[iRow][bin] = 0;
 863           if (fAllBins[iRow][bin]) fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
 864         }
 865       }
 866     }
 867   }
 868
 869   if (AliTPCReconstructor::StreamLevel()>5) {
 870     for (Int_t iRow = 0; iRow < nRows; iRow++) {
 871       Int_t maxPad = fParam->GetNPads(fSector,iRow);
 872
 873       for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
 874         for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
 875           Int_t bin = iPad*fMaxTime+iTimeBin;
 876           Float_t signal = fAllBins[iRow][bin];
 877           if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
 878             Double_t x[]={iRow,iPad-3,iTimeBin-3};
 879             Int_t i[]={fSector};
 880             AliTPCTransform trafo;
 881             trafo.Transform(x,i,0,1);
 882             Double_t gx[3]={x[0],x[1],x[2]};
 883             trafo.RotatedGlobal2Global(fSector,gx);
 884         //        fAllSigBins[iRow][fAllNSigBins[iRow]++]
 885             Int_t rowsigBins = fAllNSigBins[iRow];
 886             Int_t first=fAllSigBins[iRow][0];
 887             Int_t last= 0;
 888         //        if (rowsigBins>0) fAllSigBins[iRow][fAllNSigBins[iRow]-1];
 889
 890             if (AliTPCReconstructor::StreamLevel()>5) {
 891               (*fDebugStreamer)<<"Digits"<<
 892                 "sec="<<fSector<<
 893                 "row="<<iRow<<
 894                 "pad="<<iPad<<
 895                 "time="<<iTimeBin<<
 896                 "sig="<<signal<<
 897                 "x="<<x[0]<<
 898                 "y="<<x[1]<<
 899                 "z="<<x[2]<<
 900                 "gx="<<gx[0]<<
 901                 "gy="<<gx[1]<<
 902                 "gz="<<gx[2]<<
 903     //
 904                 "rowsigBins="<<rowsigBins<<
 905                 "first="<<first<<
 906                 "last="<<last<<
 907                 "\n";
 908             }
 909           }
 910         }
 911       }
 912     }
 913   }
 914
 915     // Now loop over rows and find clusters
 916   for (fRow = 0; fRow < nRows; fRow++) {
 917     fRowCl->SetID(fParam->GetIndex(fSector, fRow));
 918     if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
 919
 920     fRx = fParam->GetPadRowRadii(fSector, fRow);
 921     fPadLength = fParam->GetPadPitchLength(fSector, fRow);
 922     fPadWidth  = fParam->GetPadPitchWidth();
 923     fMaxPad = fParam->GetNPads(fSector,fRow);
 924     fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 925
 926     fBins = fAllBins[fRow];
 927     fSigBins = fAllSigBins[fRow];
 928     fNSigBins = fAllNSigBins[fRow];
 929
 930     FindClusters(noiseROC);
 931
 932     FillRow();
 933     if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear("C");
 934     fNclusters += fNcluster;
 935
 936   } // End of loop to find clusters
 937 }
 938
 939
 940 void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
 941 {
 942 //-----------------------------------------------------------------
 943 // This is a cluster finder for the TPC raw data.
 944 // The method assumes NO ordering of the altro channels.
 945 // The pedestal subtraction can be switched on and off
 946 // using an option of the TPC reconstructor
 947 //-----------------------------------------------------------------
 948   fRecoParam = AliTPCReconstructor::GetRecoParam();
 949   if (!fRecoParam){
 950     AliFatal("Can not get the reconstruction parameters");
 951   }
 952   if(AliTPCReconstructor::StreamLevel()>5) {
 953     AliInfo("Parameter Dumps");
 954     fParam->Dump();
 955     fRecoParam->Dump();
 956   }
 957   fRowDig = NULL;
 958
 959   //-----------------------------------------------------------------
 960   // Use HLT clusters
 961   //-----------------------------------------------------------------
 962   if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
 963     AliInfo("Using HLT clusters for TPC off-line reconstruction");
 964     fZWidth = fParam->GetZWidth();
 965     Int_t iResult = ReadHLTClusters();
 966
 967     // HLT clusters present
 968     if (iResult >= 0 && fNclusters > 0)
 969       return;
 970
 971     // HLT clusters not present
 972     if (iResult < 0 || fNclusters == 0) {
 973       if (fUseHLTClusters == 3) {
 974         AliError("No HLT clusters present, but requested.");
 975         return;
 976       }
 977       else {
 978         AliInfo("Now trying to read TPC RAW");
 979       }
 980     }
 981     // Some other problem during cluster reading
 982     else {
 983       AliWarning("Some problem while unpacking of HLT clusters.");
 984       return;
 985     }
 986   } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
 987
 988   //-----------------------------------------------------------------
 989   // Run TPC off-line clusterer
 990   //-----------------------------------------------------------------
 991   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
 992   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
 993   //
 994   AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
 995   fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
 996   if (fEventHeader){
 997     fTimeStamp = fEventHeader->Get("Timestamp");
 998     fEventType = fEventHeader->Get("Type");
 999     AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1000     transform->SetCurrentTimeStamp(fTimeStamp);
1001     transform->SetCurrentRun(rawReader->GetRunNumber());
1002   }
1003
1004   // creaate one TClonesArray for all clusters
1005   if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
1006   // reset counter
1007   fNclusters  = 0;
1008
1009   fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
1010 //   const Int_t kNIS = fParam->GetNInnerSector();
1011 //   const Int_t kNOS = fParam->GetNOuterSector();
1012 //   const Int_t kNS = kNIS + kNOS;
1013   fZWidth = fParam->GetZWidth();
1014   Int_t zeroSup = fParam->GetZeroSup();
1015   //
1016   // Clean-up
1017   //
1018   AliTPCROC * roc = AliTPCROC::Instance();
1019   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1020   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
1021   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1022     //
1023     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
1024     memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1025     fAllNSigBins[iRow]=0;
1026   }
1027
1028   Int_t prevSector=-1;
1029   rawReader->Reset();
1030   Int_t digCounter=0;
1031   //
1032   // Loop over DDLs
1033   //
1034   const Int_t kNIS = fParam->GetNInnerSector();
1035   const Int_t kNOS = fParam->GetNOuterSector();
1036   const Int_t kNS = kNIS + kNOS;
1037
1038   for(fSector = 0; fSector < kNS; fSector++) {
1039
1040     Int_t nRows = 0;
1041     Int_t nDDLs = 0, indexDDL = 0;
1042     if (fSector < kNIS) {
1043       nRows = fParam->GetNRowLow();
1044       fSign = (fSector < kNIS/2) ? 1 : -1;
1045       nDDLs = 2;
1046       indexDDL = fSector * 2;
1047     }
1048     else {
1049       nRows = fParam->GetNRowUp();
1050       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1051       nDDLs = 4;
1052       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1053     }
1054
1055     // load the raw data for corresponding DDLs
1056     rawReader->Reset();
1057     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1058
1059   while (input.NextDDL()){
1060     if (input.GetSector() != fSector)
1061       AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1062
1063     //Int_t nRows = fParam->GetNRow(fSector);
1064
1065     AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
1066     // Begin loop over altro data
1067     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1068     Float_t gain =1;
1069
1070     //loop over pads
1071     while ( input.NextChannel() ) {
1072       Int_t iRow = input.GetRow();
1073       if (iRow < 0){
1074         continue;
1075       }
1076       if (iRow >= nRows){
1077         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
1078                       iRow, 0, nRows -1));
1079         continue;
1080       }
1081       //pad
1082       Int_t iPad = input.GetPad();
1083       if (iPad < 0 || iPad >= nPadsMax) {
1084         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1085                       iPad, 0, nPadsMax-1));
1086         continue;
1087       }
1088       gain    = gainROC->GetValue(iRow,iPad);
1089       iPad+=3;
1090
1091       //loop over bunches
1092       while ( input.NextBunch() ){
1093         Int_t  startTbin    = (Int_t)input.GetStartTimeBin();
1094         Int_t  bunchlength  = (Int_t)input.GetBunchLength();
1095         const UShort_t *sig = input.GetSignals();
1096         for (Int_t iTime = 0; iTime<bunchlength; iTime++){
1097           Int_t iTimeBin=startTbin-iTime;
1098           if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1099             continue;
1100             AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1101                           iTimeBin, 0, iTimeBin -1));
1102           }
1103           iTimeBin+=3;
1104           //signal
1105           Float_t signal=(Float_t)sig[iTime];
1106           if (!calcPedestal && signal <= zeroSup) continue;
1107
1108           if (!calcPedestal) {
1109             Int_t bin = iPad*fMaxTime+iTimeBin;
1110             if (gain>0){
1111               fAllBins[iRow][bin] = signal/gain;
1112             }else{
1113               fAllBins[iRow][bin] =0;
1114             }
1115             fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
1116           }else{
1117             fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
1118           }
1119           fAllBins[iRow][iPad*fMaxTime+0]+=1.;  // pad with signal
1120
1121           // Temporary
1122           digCounter++;
1123         }// end loop signals in bunch
1124       }// end loop bunches
1125     } // end loop pads
1126     //
1127     //
1128     //
1129     //
1130     // Now loop over rows and perform pedestal subtraction
1131     if (digCounter==0) continue;
1132   } // End of loop over sectors
1133   //process last sector
1134   if ( digCounter>0 ){
1135     ProcessSectorData();
1136     for (Int_t iRow = 0; iRow < fParam->GetNRow(fSector); iRow++) {
1137       Int_t maxPad = fParam->GetNPads(fSector,iRow);
1138       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
1139       memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1140       fAllNSigBins[iRow] = 0;
1141     }
1142     prevSector=fSector;
1143     digCounter=0;
1144   }
1145   }
1146
1147   if (rawReader->GetEventId() && fOutput ){
1148     Info("Digits2Clusters", "File  %s Event\t%u\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
1149   }
1150
1151   if(rawReader->GetEventId()) {
1152     Info("Digits2Clusters", "Event\t%u\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1153   }
1154
1155   if(fBClonesArray) {
1156     //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1157   }
1158
1159   if (fUseHLTClusters == 2 && fNclusters == 0) {
1160     AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters.");
1161
1162     fZWidth = fParam->GetZWidth();
1163     ReadHLTClusters();
1164   }
1165 }
1166
1167
1168
1169
1170
1171 void AliTPCclustererMI::Digits2ClustersOld
1172 (AliRawReader* rawReader)
1173 {
1174 //-----------------------------------------------------------------
1175 // This is a cluster finder for the TPC raw data.
1176 // The method assumes NO ordering of the altro channels.
1177 // The pedestal subtraction can be switched on and off
1178 // using an option of the TPC reconstructor
1179 //-----------------------------------------------------------------
1180   fRecoParam = AliTPCReconstructor::GetRecoParam();
1181   if (!fRecoParam){
1182     AliFatal("Can not get the reconstruction parameters");
1183   }
1184   if(AliTPCReconstructor::StreamLevel()>5) {
1185     AliInfo("Parameter Dumps");
1186     fParam->Dump();
1187     fRecoParam->Dump();
1188   }
1189   fRowDig = NULL;
1190
1191   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
1192   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
1193   //
1194   AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
1195   fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1196   if (fEventHeader){
1197     fTimeStamp = fEventHeader->Get("Timestamp");
1198     fEventType = fEventHeader->Get("Type");
1199   }
1200
1201   // creaate one TClonesArray for all clusters
1202   if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
1203   // reset counter
1204   fNclusters  = 0;
1205
1206   fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
1207   const Int_t kNIS = fParam->GetNInnerSector();
1208   const Int_t kNOS = fParam->GetNOuterSector();
1209   const Int_t kNS = kNIS + kNOS;
1210   fZWidth = fParam->GetZWidth();
1211   Int_t zeroSup = fParam->GetZeroSup();
1212   //
1213   // Clean-up
1214   //
1215
1216   AliTPCROC * roc = AliTPCROC::Instance();
1217   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1218   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
1219   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1220     //
1221     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
1222     memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1223     fAllNSigBins[iRow]=0;
1224   }
1225   //
1226   // Loop over sectors
1227   //
1228   for(fSector = 0; fSector < kNS; fSector++) {
1229
1230     Int_t nRows = 0;
1231     Int_t nDDLs = 0, indexDDL = 0;
1232     if (fSector < kNIS) {
1233       nRows = fParam->GetNRowLow();
1234       fSign = (fSector < kNIS/2) ? 1 : -1;
1235       nDDLs = 2;
1236       indexDDL = fSector * 2;
1237     }
1238     else {
1239       nRows = fParam->GetNRowUp();
1240       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1241       nDDLs = 4;
1242       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1243     }
1244
1245     // load the raw data for corresponding DDLs
1246     rawReader->Reset();
1247     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1248
1249     // select only good sector
1250     if (!input.Next()) continue;
1251     if(input.GetSector() != fSector) continue;
1252
1253     AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
1254
1255     for (Int_t iRow = 0; iRow < nRows; iRow++) {
1256       Int_t maxPad;
1257       if (fSector < kNIS)
1258         maxPad = fParam->GetNPadsLow(iRow);
1259       else
1260         maxPad = fParam->GetNPadsUp(iRow);
1261
1262       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
1263       memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1264       fAllNSigBins[iRow] = 0;
1265     }
1266
1267     Int_t digCounter=0;
1268     // Begin loop over altro data
1269     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1270     Float_t gain =1;
1271     Int_t lastPad=-1;
1272
1273     input.Reset();
1274     while (input.Next()) {
1275       if (input.GetSector() != fSector)
1276         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1277
1278
1279       Int_t iRow = input.GetRow();
1280       if (iRow < 0){
1281         continue;
1282       }
1283
1284       if (iRow < 0 || iRow >= nRows){
1285         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
1286                       iRow, 0, nRows -1));
1287         continue;
1288       }
1289       //pad
1290       Int_t iPad = input.GetPad();
1291       if (iPad < 0 || iPad >= nPadsMax) {
1292         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1293                       iPad, 0, nPadsMax-1));
1294         continue;
1295       }
1296       if (iPad!=lastPad){
1297         gain    = gainROC->GetValue(iRow,iPad);
1298         lastPad = iPad;
1299       }
1300       iPad+=3;
1301       //time
1302       Int_t iTimeBin = input.GetTime();
1303       if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1304         continue;
1305         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1306                       iTimeBin, 0, iTimeBin -1));
1307       }
1308       iTimeBin+=3;
1309
1310       //signal
1311       Float_t signal = input.GetSignal();
1312       if (!calcPedestal && signal <= zeroSup) continue;
1313
1314       if (!calcPedestal) {
1315         Int_t bin = iPad*fMaxTime+iTimeBin;
1316         if (gain>0){
1317           fAllBins[iRow][bin] = signal/gain;
1318         }else{
1319           fAllBins[iRow][bin] =0;
1320         }
1321         fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
1322       }else{
1323         fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
1324       }
1325       fAllBins[iRow][iPad*fMaxTime+0]+=1.;  // pad with signal
1326
1327       // Temporary
1328       digCounter++;
1329     } // End of the loop over altro data
1330     //
1331     //
1332     //
1333     //
1334     // Now loop over rows and perform pedestal subtraction
1335     if (digCounter==0) continue;
1336     ProcessSectorData();
1337   } // End of loop over sectors
1338
1339   if (rawReader->GetEventId() && fOutput ){
1340     Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
1341   }
1342
1343   if(rawReader->GetEventId()) {
1344     Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1345   }
1346
1347   if(fBClonesArray) {
1348     //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1349   }
1350 }
1351
1352 void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
1353 {
1354
1355   //
1356   // add virtual charge at the edge
1357   //
1358   Double_t kMaxDumpSize = 500000;
1359   if (!fOutput) {
1360     fBDumpSignal =kFALSE;
1361   }else{
1362     if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE;   //dump signal flag
1363   }
1364
1365   fNcluster=0;
1366   fLoop=1;
1367   Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
1368   Float_t minMaxCutAbs       = fRecoParam->GetMinMaxCutAbs();
1369   Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1370   Float_t minUpDownCutAbs    = fRecoParam->GetMinUpDownCutAbs();
1371   Float_t minMaxCutSigma       = fRecoParam->GetMinMaxCutSigma();
1372   Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1373   Float_t minUpDownCutSigma    = fRecoParam->GetMinUpDownCutSigma();
1374   Int_t   useOnePadCluster     = fRecoParam->GetUseOnePadCluster();
1375   for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1376     Int_t i = fSigBins[iSig];
1377     if (i%fMaxTime<=crtime) continue;
1378     Float_t *b = &fBins[i];
1379     //absolute custs
1380     if (b[0]<minMaxCutAbs) continue;   //threshold for maxima
1381     //
1382     if (useOnePadCluster==0){
1383       if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue;  // cut on isolated clusters
1384       if (b[-1]+b[1]<=0) continue;               // cut on isolated clusters
1385       if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1386     }
1387     //
1388     if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue;   //threshold for up down  (TRF)
1389     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue;   //threshold for left right (PRF)
1390     if (!IsMaximum(*b,fMaxTime,b)) continue;
1391     //
1392     Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
1393     if (noise>fRecoParam->GetMaxNoise()) continue;
1394     // sigma cuts
1395     if (b[0]<minMaxCutSigma*noise) continue;   //threshold form maxima
1396     if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue;   //threshold for up town TRF
1397     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue;   //threshold for left right (PRF)
1398
1399     AliTPCclusterMI c;   // default cosntruction  without info
1400     Int_t dummy=0;
1401     MakeCluster(i, fMaxTime, fBins, dummy,c);
1402
1403     //}
1404   }
1405 }
1406
1407 Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
1408   // -- Depricated --
1409   // Currently hack to filter digital noise (15.06.2008)
1410   // To be parameterized in the AliTPCrecoParam
1411   // More inteligent way  to be used in future
1412   // Acces to the proper pedestal file needed
1413   //
1414   if (cl->GetMax()<400) return kTRUE;
1415   Double_t ratio = cl->GetQ()/cl->GetMax();
1416   if (cl->GetMax()>700){
1417     if ((ratio - int(ratio)>0.8)) return kFALSE;
1418   }
1419   if ((ratio - int(ratio)<0.95)) return kTRUE;
1420   return kFALSE;
1421 }
1422
1423
1424 Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
1425   //
1426   // process signal on given pad - + streaming of additional information in special mode
1427   //
1428   // id[0] - sector
1429   // id[1] - row
1430   // id[2] - pad
1431
1432   //
1433   // ESTIMATE pedestal and the noise
1434   //
1435   const Int_t kPedMax = 100;
1436   Float_t  max    =  0;
1437   Float_t  maxPos =  0;
1438   Int_t    median =  -1;
1439   Int_t    count0 =  0;
1440   Int_t    count1 =  0;
1441   Float_t  rmsCalib   = rmsEvent;       // backup initial value ( from calib)
1442   Float_t  pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1443   Int_t    firstBin = fRecoParam->GetFirstBin();
1444   //
1445   UShort_t histo[kPedMax];
1446   //memset(histo,0,kPedMax*sizeof(UShort_t));
1447   for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
1448   for (Int_t i=0; i<fMaxTime; i++){
1449     if (signal[i]<=0) continue;
1450     if (signal[i]>max && i>firstBin) {
1451       max = signal[i];
1452       maxPos = i;
1453     }
1454     if (signal[i]>kPedMax-1) continue;
1455     histo[int(signal[i]+0.5)]++;
1456     count0++;
1457   }
1458   //
1459   for (Int_t i=1; i<kPedMax; i++){
1460     if (count1<count0*0.5) median=i;
1461     count1+=histo[i];
1462   }
1463   // truncated mean
1464   //
1465   Float_t count10=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
1466   Float_t count06=histo[median] ,mean06=histo[median]*median,  rms06=histo[median]*median*median ;
1467   Float_t count09=histo[median] ,mean09=histo[median]*median,  rms09=histo[median]*median*median ;
1468   //
1469   for (Int_t idelta=1; idelta<10; idelta++){
1470     if (median-idelta<=0) continue;
1471     if (median+idelta>kPedMax) continue;
1472     if (count06<0.6*count1){
1473       count06+=histo[median-idelta];
1474       mean06 +=histo[median-idelta]*(median-idelta);
1475       rms06  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1476       count06+=histo[median+idelta];
1477       mean06 +=histo[median+idelta]*(median+idelta);
1478       rms06  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1479     }
1480     if (count09<0.9*count1){
1481       count09+=histo[median-idelta];
1482       mean09 +=histo[median-idelta]*(median-idelta);
1483       rms09  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1484       count09+=histo[median+idelta];
1485       mean09 +=histo[median+idelta]*(median+idelta);
1486       rms09  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1487     }
1488     if (count10<0.95*count1){
1489       count10+=histo[median-idelta];
1490       mean +=histo[median-idelta]*(median-idelta);
1491       rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1492       count10+=histo[median+idelta];
1493       mean +=histo[median+idelta]*(median+idelta);
1494       rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1495     }
1496   }
1497   if (count10) {
1498     mean  /=count10;
1499     rms    = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1500   }
1501   if (count06) {
1502     mean06/=count06;
1503     rms06    = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1504   }
1505   if (count09) {
1506     mean09/=count09;
1507     rms09    = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1508   }
1509   rmsEvent = rms09;
1510   //
1511   pedestalEvent = median;
1512   if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1513   //
1514   UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1515   //
1516   // Dump mean signal info
1517   //
1518     if (AliTPCReconstructor::StreamLevel()>0) {
1519     (*fDebugStreamer)<<"Signal"<<
1520     "TimeStamp="<<fTimeStamp<<
1521     "EventType="<<fEventType<<
1522     "Sector="<<uid[0]<<
1523     "Row="<<uid[1]<<
1524     "Pad="<<uid[2]<<
1525     "Max="<<max<<
1526     "MaxPos="<<maxPos<<
1527     //
1528     "Median="<<median<<
1529     "Mean="<<mean<<
1530     "RMS="<<rms<<
1531     "Mean06="<<mean06<<
1532     "RMS06="<<rms06<<
1533     "Mean09="<<mean09<<
1534     "RMS09="<<rms09<<
1535     "RMSCalib="<<rmsCalib<<
1536     "PedCalib="<<pedestalCalib<<
1537     "\n";
1538     }
1539   //
1540   // fill pedestal histogram
1541   //
1542   //
1543   //
1544   //
1545   Float_t kMin =fRecoParam->GetDumpAmplitudeMin();   // minimal signal to be dumped
1546   Float_t *dsignal = new Float_t[nchannels];
1547   Float_t *dtime   = new Float_t[nchannels];
1548   for (Int_t i=0; i<nchannels; i++){
1549     dtime[i] = i;
1550     dsignal[i] = signal[i];
1551   }
1552
1553   //
1554   // Big signals dumping
1555   //
1556   if (AliTPCReconstructor::StreamLevel()>0) {
1557   if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin())
1558     (*fDebugStreamer)<<"SignalB"<<     // pads with signal
1559       "TimeStamp="<<fTimeStamp<<
1560       "EventType="<<fEventType<<
1561       "Sector="<<uid[0]<<
1562       "Row="<<uid[1]<<
1563       "Pad="<<uid[2]<<
1564       //      "Graph="<<graph<<
1565       "Max="<<max<<
1566       "MaxPos="<<maxPos<<
1567       //
1568       "Median="<<median<<
1569       "Mean="<<mean<<
1570       "RMS="<<rms<<
1571       "Mean06="<<mean06<<
1572       "RMS06="<<rms06<<
1573       "Mean09="<<mean09<<
1574       "RMS09="<<rms09<<
1575       "\n";
1576   }
1577
1578   delete [] dsignal;
1579   delete [] dtime;
1580   if (rms06>fRecoParam->GetMaxNoise()) {
1581     pedestalEvent+=1024.;
1582     return 1024+median; // sign noisy channel in debug mode
1583   }
1584   return median;
1585 }
1586
1587 Int_t AliTPCclustererMI::ReadHLTClusters()
1588 {
1589   //
1590   // read HLT clusters instead of off line custers,
1591   // used in Digits2Clusters
1592   //
1593
1594   if (!fHLTClusterAccess) {
1595   TClass* pCl=NULL;
1596   ROOT::NewFunc_t pNewFunc=NULL;
1597   do {
1598     pCl=TClass::GetClass("AliHLTTPCClusterAccessHLTOUT");
1599   } while (!pCl && gSystem->Load("libAliHLTTPC.so")==0);
1600   if (!pCl || (pNewFunc=pCl->GetNew())==NULL) {
1601     AliError("can not load class description of AliHLTTPCClusterAccessHLTOUT, aborting ...");
1602     return -1;
1603   }
1604
1605   void* p=(*pNewFunc)(NULL);
1606   if (!p) {
1607     AliError("unable to create instance of AliHLTTPCClusterAccessHLTOUT");
1608     return -2;
1609   }
1610   fHLTClusterAccess=reinterpret_cast<TObject*>(p);
1611   }
1612
1613   TObject* pClusterAccess=fHLTClusterAccess;
1614
1615   const Int_t kNIS = fParam->GetNInnerSector();
1616   const Int_t kNOS = fParam->GetNOuterSector();
1617   const Int_t kNS = kNIS + kNOS;
1618   fNclusters  = 0;
1619
1620   // make sure that all clusters from the previous event are cleared
1621   pClusterAccess->Clear("event");
1622   for(fSector = 0; fSector < kNS; fSector++) {
1623
1624     Int_t iResult = 1;
1625     TString param("sector="); param+=fSector;
1626     // prepare for next sector
1627     pClusterAccess->Clear("sector");
1628     pClusterAccess->Execute("read", param, &iResult);
1629     if (iResult < 0) {
1630       return iResult;
1631       AliError("HLT Clusters can not be found");
1632     }
1633
1634     TObject* pObj=pClusterAccess->FindObject("clusterarray");
1635     if (pObj==NULL) {
1636       AliError("HLT clusters requested, but not cluster array not present");
1637       return -4;
1638     }
1639
1640     TObjArray* clusterArray=dynamic_cast<TClonesArray*>(pObj);
1641     if (!clusterArray) {
1642       AliError("HLT cluster array is not of class type TClonesArray");
1643       return -5;
1644     }
1645
1646     AliDebug(4,Form("Reading %d clusters from HLT for sector %d", clusterArray->GetEntriesFast(), fSector));
1647
1648     Int_t nClusterSector=0;
1649     Int_t nRows=fParam->GetNRow(fSector);
1650
1651     for (fRow = 0; fRow < nRows; fRow++) {
1652       fRowCl->SetID(fParam->GetIndex(fSector, fRow));
1653       if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
1654       fNcluster=0; // reset clusters per row
1655
1656       fRx = fParam->GetPadRowRadii(fSector, fRow);
1657       fPadLength = fParam->GetPadPitchLength(fSector, fRow);
1658       fPadWidth  = fParam->GetPadPitchWidth();
1659       fMaxPad = fParam->GetNPads(fSector,fRow);
1660       fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
1661
1662       fBins = fAllBins[fRow];
1663       fSigBins = fAllSigBins[fRow];
1664       fNSigBins = fAllNSigBins[fRow];
1665
1666       for (Int_t i=0; i<clusterArray->GetEntriesFast(); i++) {
1667         if (!clusterArray->At(i))
1668           continue;
1669
1670         AliTPCclusterMI* cluster=dynamic_cast<AliTPCclusterMI*>(clusterArray->At(i));
1671         if (!cluster) continue;
1672         if (cluster->GetRow()!=fRow) continue;
1673         nClusterSector++;
1674         AddCluster(*cluster, NULL, 0);
1675       }
1676
1677       FillRow();
1678       fRowCl->GetArray()->Clear("c");
1679
1680     } // for (fRow = 0; fRow < nRows; fRow++) {
1681     if (nClusterSector!=clusterArray->GetEntriesFast()) {
1682       AliError(Form("Failed to read %d out of %d HLT clusters",
1683                     clusterArray->GetEntriesFast()-nClusterSector,
1684                     clusterArray->GetEntriesFast()));
1685     }
1686     fNclusters+=nClusterSector;
1687   } // for(fSector = 0; fSector < kNS; fSector++) {
1688
1689   pClusterAccess->Clear("event");
1690
1691   Info("Digits2Clusters", "Number of converted HLT clusters : %d", fNclusters);
1692
1693   return 0;
1694 }