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