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   AliInfo(Form("Usage of HLT clusters in TPC reconstruction : %d",fUseHLTClusters));
 731
 732   if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
 733     AliInfo("Using HLT clusters for TPC off-line reconstruction");
 734     fZWidth = fParam->GetZWidth();
 735     Int_t iResult = ReadHLTClusters();
 736
 737     // HLT clusters present
 738     if (!iResult)
 739       return;
 740     // HLT clusters not present
 741     else if(iResult == -1) {
 742       if (fUseHLTClusters == 3) {
 743         AliError("No HLT clusters present, but requiered.");
 744         return;
 745       }
 746     }
 747     // Some other problem during cluster reading
 748     else {
 749       AliWarning("Some problem while unpacking of HLT clusters.");
 750       return;
 751     }
 752   } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
 753
 754   //-----------------------------------------------------------------
 755   // Run TPC off-line clusterer
 756   //-----------------------------------------------------------------
 757   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
 758   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 759   AliSimDigits digarr, *dummy=&digarr;
 760   fRowDig = dummy;
 761   fInput->GetBranch("Segment")->SetAddress(&dummy);
 762   Stat_t nentries = fInput->GetEntries();
 763
 764   fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3   after
 765
 766   Int_t nclusters  = 0;
 767
 768   for (Int_t n=0; n<nentries; n++) {
 769     fInput->GetEvent(n);
 770     if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
 771       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
 772       continue;
 773     }
 774     Int_t row = fRow;
 775     AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector);  // pad gains per given sector
 776     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector); // noise per given sector
 777     //
 778
 779     fRowCl->SetID(digarr.GetID());
 780     if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
 781     fRx=fParam->GetPadRowRadii(fSector,row);
 782
 783
 784     const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
 785     fZWidth = fParam->GetZWidth();
 786     if (fSector < kNIS) {
 787       fMaxPad = fParam->GetNPadsLow(row);
 788       fSign = (fSector < kNIS/2) ? 1 : -1;
 789       fPadLength = fParam->GetPadPitchLength(fSector,row);
 790       fPadWidth = fParam->GetPadPitchWidth();
 791     } else {
 792       fMaxPad = fParam->GetNPadsUp(row);
 793       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 794       fPadLength = fParam->GetPadPitchLength(fSector,row);
 795       fPadWidth  = fParam->GetPadPitchWidth();
 796     }
 797
 798
 799     fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 800     fBins    =new Float_t[fMaxBin];
 801     fSigBins =new Int_t[fMaxBin];
 802     fNSigBins = 0;
 803     memset(fBins,0,sizeof(Float_t)*fMaxBin);
 804
 805     if (digarr.First()) //MI change
 806       do {
 807         Float_t dig=digarr.CurrentDigit();
 808         if (dig<=fParam->GetZeroSup()) continue;
 809         Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
 810         Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
 811         Int_t bin = i*fMaxTime+j;
 812         if (gain>0){
 813           fBins[bin]=dig/gain;
 814         }else{
 815           fBins[bin]=0;
 816         }
 817         fSigBins[fNSigBins++]=bin;
 818       } while (digarr.Next());
 819     digarr.ExpandTrackBuffer();
 820
 821     FindClusters(noiseROC);
 822     FillRow();
 823     fRowCl->GetArray()->Clear("C");
 824     nclusters+=fNcluster;
 825
 826     delete[] fBins;
 827     delete[] fSigBins;
 828   }
 829
 830   Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
 831 }
 832
 833 void AliTPCclustererMI::ProcessSectorData(){
 834   //
 835   // Process the data for the current sector
 836   //
 837
 838   AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
 839   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 840   AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector);  // pedestal per given sector
 841   AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector);  // noise per given sector
 842   //check the presence of the calibration
 843   if (!noiseROC ||!pedestalROC ) {
 844     AliError(Form("Missing calibration per sector\t%d\n",fSector));
 845     return;
 846   }
 847   Int_t  nRows=fParam->GetNRow(fSector);
 848   Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
 849   Int_t zeroSup = fParam->GetZeroSup();
 850   //    if (calcPedestal) {
 851   if (kFALSE ) {
 852     for (Int_t iRow = 0; iRow < nRows; iRow++) {
 853       Int_t maxPad = fParam->GetNPads(fSector, iRow);
 854
 855       for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
 856     //
 857     // Temporary fix for data production - !!!! MARIAN
 858     // The noise calibration should take mean and RMS - currently the Gaussian fit used
 859     // In case of double peak  - the pad should be rejected
 860     //
 861     // Line mean - if more than given digits over threshold - make a noise calculation
 862     // and pedestal substration
 863         if (!calcPedestal && fAllBins[iRow][iPad*fMaxTime+0]<50) continue;
 864     //
 865         if (fAllBins[iRow][iPad*fMaxTime+0] <1 ) continue;  // no data
 866         Float_t *p = &fAllBins[iRow][iPad*fMaxTime+3];
 867     //Float_t pedestal = TMath::Median(fMaxTime, p);
 868         Int_t id[3] = {fSector, iRow, iPad-3};
 869     // calib values
 870         Double_t rmsCalib=  noiseROC->GetValue(iRow,iPad-3);
 871         Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
 872         Double_t rmsEvent       = rmsCalib;
 873         Double_t pedestalEvent  = pedestalCalib;
 874         ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
 875         if (rmsEvent<rmsCalib) rmsEvent = rmsCalib;   // take worst scenario
 876         if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
 877
 878     //
 879         for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
 880           Int_t bin = iPad*fMaxTime+iTimeBin;
 881           fAllBins[iRow][bin] -= pedestalEvent;
 882           if (iTimeBin < fRecoParam->GetFirstBin())
 883             fAllBins[iRow][bin] = 0;
 884           if (iTimeBin > fRecoParam->GetLastBin())
 885             fAllBins[iRow][bin] = 0;
 886           if (fAllBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
 887             fAllBins[iRow][bin] = 0;
 888           if (fAllBins[iRow][bin] < 3.0*rmsEvent)   // 3 sigma cut on RMS
 889             fAllBins[iRow][bin] = 0;
 890           if (fAllBins[iRow][bin]) fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
 891         }
 892       }
 893     }
 894   }
 895
 896   if (AliTPCReconstructor::StreamLevel()>5) {
 897     for (Int_t iRow = 0; iRow < nRows; iRow++) {
 898       Int_t maxPad = fParam->GetNPads(fSector,iRow);
 899
 900       for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
 901         for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
 902           Int_t bin = iPad*fMaxTime+iTimeBin;
 903           Float_t signal = fAllBins[iRow][bin];
 904           if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
 905             Double_t x[]={iRow,iPad-3,iTimeBin-3};
 906             Int_t i[]={fSector};
 907             AliTPCTransform trafo;
 908             trafo.Transform(x,i,0,1);
 909             Double_t gx[3]={x[0],x[1],x[2]};
 910             trafo.RotatedGlobal2Global(fSector,gx);
 911         //        fAllSigBins[iRow][fAllNSigBins[iRow]++]
 912             Int_t rowsigBins = fAllNSigBins[iRow];
 913             Int_t first=fAllSigBins[iRow][0];
 914             Int_t last= 0;
 915         //        if (rowsigBins>0) fAllSigBins[iRow][fAllNSigBins[iRow]-1];
 916
 917             if (AliTPCReconstructor::StreamLevel()>5) {
 918               (*fDebugStreamer)<<"Digits"<<
 919                 "sec="<<fSector<<
 920                 "row="<<iRow<<
 921                 "pad="<<iPad<<
 922                 "time="<<iTimeBin<<
 923                 "sig="<<signal<<
 924                 "x="<<x[0]<<
 925                 "y="<<x[1]<<
 926                 "z="<<x[2]<<
 927                 "gx="<<gx[0]<<
 928                 "gy="<<gx[1]<<
 929                 "gz="<<gx[2]<<
 930     //
 931                 "rowsigBins="<<rowsigBins<<
 932                 "first="<<first<<
 933                 "last="<<last<<
 934                 "\n";
 935             }
 936           }
 937         }
 938       }
 939     }
 940   }
 941
 942     // Now loop over rows and find clusters
 943   for (fRow = 0; fRow < nRows; fRow++) {
 944     fRowCl->SetID(fParam->GetIndex(fSector, fRow));
 945     if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
 946
 947     fRx = fParam->GetPadRowRadii(fSector, fRow);
 948     fPadLength = fParam->GetPadPitchLength(fSector, fRow);
 949     fPadWidth  = fParam->GetPadPitchWidth();
 950     fMaxPad = fParam->GetNPads(fSector,fRow);
 951     fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 952
 953     fBins = fAllBins[fRow];
 954     fSigBins = fAllSigBins[fRow];
 955     fNSigBins = fAllNSigBins[fRow];
 956
 957     FindClusters(noiseROC);
 958
 959     FillRow();
 960     if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear("C");
 961     fNclusters += fNcluster;
 962
 963   } // End of loop to find clusters
 964 }
 965
 966
 967 void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
 968 {
 969 //-----------------------------------------------------------------
 970 // This is a cluster finder for the TPC raw data.
 971 // The method assumes NO ordering of the altro channels.
 972 // The pedestal subtraction can be switched on and off
 973 // using an option of the TPC reconstructor
 974 //-----------------------------------------------------------------
 975   fRecoParam = AliTPCReconstructor::GetRecoParam();
 976   if (!fRecoParam){
 977     AliFatal("Can not get the reconstruction parameters");
 978   }
 979   if(AliTPCReconstructor::StreamLevel()>5) {
 980     AliInfo("Parameter Dumps");
 981     fParam->Dump();
 982     fRecoParam->Dump();
 983   }
 984   fRowDig = NULL;
 985
 986   //-----------------------------------------------------------------
 987   // Use HLT clusters
 988   //-----------------------------------------------------------------
 989   fUseHLTClusters = fRecoParam->GetUseHLTClusters();
 990
 991   AliInfo(Form("Usage of HLT clusters in TPC reconstruction : %d",fUseHLTClusters));
 992
 993   if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
 994     AliInfo("Using HLT clusters for TPC off-line reconstruction");
 995     fZWidth = fParam->GetZWidth();
 996     Int_t iResult = ReadHLTClusters();
 997
 998     AliError(Form("HLT result : %d",iResult));
 999
1000     // HLT clusters present
1001     if (!iResult)
1002       return;
1003     // HLT clusters not present
1004     else if(iResult == -1) {
1005       if (fUseHLTClusters == 3) {
1006         AliError("No HLT clusters present, but requiered.");
1007         return;
1008       }
1009     }
1010     // Some other problem during cluster reading
1011     else {
1012       AliWarning("Some problem while unpacking of HLT clusters.");
1013       return;
1014     }
1015   } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
1016
1017   //-----------------------------------------------------------------
1018   // Run TPC off-line clusterer
1019   //-----------------------------------------------------------------
1020   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
1021   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
1022   //
1023   AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
1024   fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1025   if (fEventHeader){
1026     fTimeStamp = fEventHeader->Get("Timestamp");
1027     fEventType = fEventHeader->Get("Type");
1028     AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1029     transform->SetCurrentTimeStamp(fTimeStamp);
1030     transform->SetCurrentRun(rawReader->GetRunNumber());
1031   }
1032
1033   // creaate one TClonesArray for all clusters
1034   if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
1035   // reset counter
1036   fNclusters  = 0;
1037
1038   fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
1039 //   const Int_t kNIS = fParam->GetNInnerSector();
1040 //   const Int_t kNOS = fParam->GetNOuterSector();
1041 //   const Int_t kNS = kNIS + kNOS;
1042   fZWidth = fParam->GetZWidth();
1043   Int_t zeroSup = fParam->GetZeroSup();
1044   //
1045   // Clean-up
1046   //
1047   AliTPCROC * roc = AliTPCROC::Instance();
1048   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1049   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
1050   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1051     //
1052     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
1053     memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1054     fAllNSigBins[iRow]=0;
1055   }
1056
1057   Int_t prevSector=-1;
1058   rawReader->Reset();
1059   Int_t digCounter=0;
1060   //
1061   // Loop over DDLs
1062   //
1063   const Int_t kNIS = fParam->GetNInnerSector();
1064   const Int_t kNOS = fParam->GetNOuterSector();
1065   const Int_t kNS = kNIS + kNOS;
1066
1067   for(fSector = 0; fSector < kNS; fSector++) {
1068
1069     Int_t nRows = 0;
1070     Int_t nDDLs = 0, indexDDL = 0;
1071     if (fSector < kNIS) {
1072       nRows = fParam->GetNRowLow();
1073       fSign = (fSector < kNIS/2) ? 1 : -1;
1074       nDDLs = 2;
1075       indexDDL = fSector * 2;
1076     }
1077     else {
1078       nRows = fParam->GetNRowUp();
1079       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1080       nDDLs = 4;
1081       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1082     }
1083
1084     // load the raw data for corresponding DDLs
1085     rawReader->Reset();
1086     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1087
1088   while (input.NextDDL()){
1089     if (input.GetSector() != fSector)
1090       AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1091
1092     //Int_t nRows = fParam->GetNRow(fSector);
1093
1094     AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
1095     // Begin loop over altro data
1096     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1097     Float_t gain =1;
1098
1099     //loop over pads
1100     while ( input.NextChannel() ) {
1101       Int_t iRow = input.GetRow();
1102       if (iRow < 0){
1103         continue;
1104       }
1105       if (iRow >= nRows){
1106         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
1107                       iRow, 0, nRows -1));
1108         continue;
1109       }
1110       //pad
1111       Int_t iPad = input.GetPad();
1112       if (iPad < 0 || iPad >= nPadsMax) {
1113         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1114                       iPad, 0, nPadsMax-1));
1115         continue;
1116       }
1117       gain    = gainROC->GetValue(iRow,iPad);
1118       iPad+=3;
1119
1120       //loop over bunches
1121       while ( input.NextBunch() ){
1122         Int_t  startTbin    = (Int_t)input.GetStartTimeBin();
1123         Int_t  bunchlength  = (Int_t)input.GetBunchLength();
1124         const UShort_t *sig = input.GetSignals();
1125         for (Int_t iTime = 0; iTime<bunchlength; iTime++){
1126           Int_t iTimeBin=startTbin-iTime;
1127           if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1128             continue;
1129             AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1130                           iTimeBin, 0, iTimeBin -1));
1131           }
1132           iTimeBin+=3;
1133           //signal
1134           Float_t signal=(Float_t)sig[iTime];
1135           if (!calcPedestal && signal <= zeroSup) continue;
1136
1137           if (!calcPedestal) {
1138             Int_t bin = iPad*fMaxTime+iTimeBin;
1139             if (gain>0){
1140               fAllBins[iRow][bin] = signal/gain;
1141             }else{
1142               fAllBins[iRow][bin] =0;
1143             }
1144             fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
1145           }else{
1146             fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
1147           }
1148           fAllBins[iRow][iPad*fMaxTime+0]+=1.;  // pad with signal
1149
1150           // Temporary
1151           digCounter++;
1152         }// end loop signals in bunch
1153       }// end loop bunches
1154     } // end loop pads
1155     //
1156     //
1157     //
1158     //
1159     // Now loop over rows and perform pedestal subtraction
1160     if (digCounter==0) continue;
1161   } // End of loop over sectors
1162   //process last sector
1163   if ( digCounter>0 ){
1164     ProcessSectorData();
1165     for (Int_t iRow = 0; iRow < fParam->GetNRow(fSector); iRow++) {
1166       Int_t maxPad = fParam->GetNPads(fSector,iRow);
1167       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
1168       memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1169       fAllNSigBins[iRow] = 0;
1170     }
1171     prevSector=fSector;
1172     digCounter=0;
1173   }
1174   }
1175
1176   if (rawReader->GetEventId() && fOutput ){
1177     Info("Digits2Clusters", "File  %s Event\t%u\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
1178   }
1179
1180   if(rawReader->GetEventId()) {
1181     Info("Digits2Clusters", "Event\t%u\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1182   }
1183
1184   if(fBClonesArray) {
1185     //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1186   }
1187 }
1188
1189
1190
1191
1192
1193 void AliTPCclustererMI::Digits2ClustersOld
1194 (AliRawReader* rawReader)
1195 {
1196 //-----------------------------------------------------------------
1197 // This is a cluster finder for the TPC raw data.
1198 // The method assumes NO ordering of the altro channels.
1199 // The pedestal subtraction can be switched on and off
1200 // using an option of the TPC reconstructor
1201 //-----------------------------------------------------------------
1202   fRecoParam = AliTPCReconstructor::GetRecoParam();
1203   if (!fRecoParam){
1204     AliFatal("Can not get the reconstruction parameters");
1205   }
1206   if(AliTPCReconstructor::StreamLevel()>5) {
1207     AliInfo("Parameter Dumps");
1208     fParam->Dump();
1209     fRecoParam->Dump();
1210   }
1211   fRowDig = NULL;
1212
1213   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
1214   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
1215   //
1216   AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
1217   fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1218   if (fEventHeader){
1219     fTimeStamp = fEventHeader->Get("Timestamp");
1220     fEventType = fEventHeader->Get("Type");
1221   }
1222
1223   // creaate one TClonesArray for all clusters
1224   if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
1225   // reset counter
1226   fNclusters  = 0;
1227
1228   fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
1229   const Int_t kNIS = fParam->GetNInnerSector();
1230   const Int_t kNOS = fParam->GetNOuterSector();
1231   const Int_t kNS = kNIS + kNOS;
1232   fZWidth = fParam->GetZWidth();
1233   Int_t zeroSup = fParam->GetZeroSup();
1234   //
1235   // Clean-up
1236   //
1237
1238   AliTPCROC * roc = AliTPCROC::Instance();
1239   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1240   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
1241   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1242     //
1243     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
1244     memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1245     fAllNSigBins[iRow]=0;
1246   }
1247   //
1248   // Loop over sectors
1249   //
1250   for(fSector = 0; fSector < kNS; fSector++) {
1251
1252     Int_t nRows = 0;
1253     Int_t nDDLs = 0, indexDDL = 0;
1254     if (fSector < kNIS) {
1255       nRows = fParam->GetNRowLow();
1256       fSign = (fSector < kNIS/2) ? 1 : -1;
1257       nDDLs = 2;
1258       indexDDL = fSector * 2;
1259     }
1260     else {
1261       nRows = fParam->GetNRowUp();
1262       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1263       nDDLs = 4;
1264       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1265     }
1266
1267     // load the raw data for corresponding DDLs
1268     rawReader->Reset();
1269     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1270
1271     // select only good sector
1272     if (!input.Next()) continue;
1273     if(input.GetSector() != fSector) continue;
1274
1275     AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
1276
1277     for (Int_t iRow = 0; iRow < nRows; iRow++) {
1278       Int_t maxPad;
1279       if (fSector < kNIS)
1280         maxPad = fParam->GetNPadsLow(iRow);
1281       else
1282         maxPad = fParam->GetNPadsUp(iRow);
1283
1284       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
1285       memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
1286       fAllNSigBins[iRow] = 0;
1287     }
1288
1289     Int_t digCounter=0;
1290     // Begin loop over altro data
1291     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1292     Float_t gain =1;
1293     Int_t lastPad=-1;
1294
1295     input.Reset();
1296     while (input.Next()) {
1297       if (input.GetSector() != fSector)
1298         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1299
1300
1301       Int_t iRow = input.GetRow();
1302       if (iRow < 0){
1303         continue;
1304       }
1305
1306       if (iRow < 0 || iRow >= nRows){
1307         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
1308                       iRow, 0, nRows -1));
1309         continue;
1310       }
1311       //pad
1312       Int_t iPad = input.GetPad();
1313       if (iPad < 0 || iPad >= nPadsMax) {
1314         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1315                       iPad, 0, nPadsMax-1));
1316         continue;
1317       }
1318       if (iPad!=lastPad){
1319         gain    = gainROC->GetValue(iRow,iPad);
1320         lastPad = iPad;
1321       }
1322       iPad+=3;
1323       //time
1324       Int_t iTimeBin = input.GetTime();
1325       if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1326         continue;
1327         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1328                       iTimeBin, 0, iTimeBin -1));
1329       }
1330       iTimeBin+=3;
1331
1332       //signal
1333       Float_t signal = input.GetSignal();
1334       if (!calcPedestal && signal <= zeroSup) continue;
1335
1336       if (!calcPedestal) {
1337         Int_t bin = iPad*fMaxTime+iTimeBin;
1338         if (gain>0){
1339           fAllBins[iRow][bin] = signal/gain;
1340         }else{
1341           fAllBins[iRow][bin] =0;
1342         }
1343         fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
1344       }else{
1345         fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
1346       }
1347       fAllBins[iRow][iPad*fMaxTime+0]+=1.;  // pad with signal
1348
1349       // Temporary
1350       digCounter++;
1351     } // End of the loop over altro data
1352     //
1353     //
1354     //
1355     //
1356     // Now loop over rows and perform pedestal subtraction
1357     if (digCounter==0) continue;
1358     ProcessSectorData();
1359   } // End of loop over sectors
1360
1361   if (rawReader->GetEventId() && fOutput ){
1362     Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
1363   }
1364
1365   if(rawReader->GetEventId()) {
1366     Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1367   }
1368
1369   if(fBClonesArray) {
1370     //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1371   }
1372 }
1373
1374 void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
1375 {
1376
1377   //
1378   // add virtual charge at the edge
1379   //
1380   Double_t kMaxDumpSize = 500000;
1381   if (!fOutput) {
1382     fBDumpSignal =kFALSE;
1383   }else{
1384     if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE;   //dump signal flag
1385   }
1386
1387   fNcluster=0;
1388   fLoop=1;
1389   Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
1390   Float_t minMaxCutAbs       = fRecoParam->GetMinMaxCutAbs();
1391   Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1392   Float_t minUpDownCutAbs    = fRecoParam->GetMinUpDownCutAbs();
1393   Float_t minMaxCutSigma       = fRecoParam->GetMinMaxCutSigma();
1394   Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1395   Float_t minUpDownCutSigma    = fRecoParam->GetMinUpDownCutSigma();
1396   Int_t   useOnePadCluster     = fRecoParam->GetUseOnePadCluster();
1397   for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1398     Int_t i = fSigBins[iSig];
1399     if (i%fMaxTime<=crtime) continue;
1400     Float_t *b = &fBins[i];
1401     //absolute custs
1402     if (b[0]<minMaxCutAbs) continue;   //threshold for maxima
1403     //
1404     if (useOnePadCluster==0){
1405       if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue;  // cut on isolated clusters
1406       if (b[-1]+b[1]<=0) continue;               // cut on isolated clusters
1407       if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1408     }
1409     //
1410     if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue;   //threshold for up down  (TRF)
1411     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue;   //threshold for left right (PRF)
1412     if (!IsMaximum(*b,fMaxTime,b)) continue;
1413     //
1414     Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
1415     if (noise>fRecoParam->GetMaxNoise()) continue;
1416     // sigma cuts
1417     if (b[0]<minMaxCutSigma*noise) continue;   //threshold form maxima
1418     if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue;   //threshold for up town TRF
1419     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue;   //threshold for left right (PRF)
1420
1421     AliTPCclusterMI c;   // default cosntruction  without info
1422     Int_t dummy=0;
1423     MakeCluster(i, fMaxTime, fBins, dummy,c);
1424
1425     //}
1426   }
1427 }
1428
1429 Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
1430   //
1431   // Currently hack to filter digital noise (15.06.2008)
1432   // To be parameterized in the AliTPCrecoParam
1433   // More inteligent way  to be used in future
1434   // Acces to the proper pedestal file needed
1435   //
1436   if (cl->GetMax()<400) return kTRUE;
1437   Double_t ratio = cl->GetQ()/cl->GetMax();
1438   if (cl->GetMax()>700){
1439     if ((ratio - int(ratio)>0.8)) return kFALSE;
1440   }
1441   if ((ratio - int(ratio)<0.95)) return kTRUE;
1442   return kFALSE;
1443 }
1444
1445
1446 Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
1447   //
1448   // process signal on given pad - + streaming of additional information in special mode
1449   //
1450   // id[0] - sector
1451   // id[1] - row
1452   // id[2] - pad
1453
1454   //
1455   // ESTIMATE pedestal and the noise
1456   //
1457   const Int_t kPedMax = 100;
1458   Float_t  max    =  0;
1459   Float_t  maxPos =  0;
1460   Int_t    median =  -1;
1461   Int_t    count0 =  0;
1462   Int_t    count1 =  0;
1463   Float_t  rmsCalib   = rmsEvent;       // backup initial value ( from calib)
1464   Float_t  pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1465   Int_t    firstBin = fRecoParam->GetFirstBin();
1466   //
1467   UShort_t histo[kPedMax];
1468   //memset(histo,0,kPedMax*sizeof(UShort_t));
1469   for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
1470   for (Int_t i=0; i<fMaxTime; i++){
1471     if (signal[i]<=0) continue;
1472     if (signal[i]>max && i>firstBin) {
1473       max = signal[i];
1474       maxPos = i;
1475     }
1476     if (signal[i]>kPedMax-1) continue;
1477     histo[int(signal[i]+0.5)]++;
1478     count0++;
1479   }
1480   //
1481   for (Int_t i=1; i<kPedMax; i++){
1482     if (count1<count0*0.5) median=i;
1483     count1+=histo[i];
1484   }
1485   // truncated mean
1486   //
1487   Float_t count10=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
1488   Float_t count06=histo[median] ,mean06=histo[median]*median,  rms06=histo[median]*median*median ;
1489   Float_t count09=histo[median] ,mean09=histo[median]*median,  rms09=histo[median]*median*median ;
1490   //
1491   for (Int_t idelta=1; idelta<10; idelta++){
1492     if (median-idelta<=0) continue;
1493     if (median+idelta>kPedMax) continue;
1494     if (count06<0.6*count1){
1495       count06+=histo[median-idelta];
1496       mean06 +=histo[median-idelta]*(median-idelta);
1497       rms06  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1498       count06+=histo[median+idelta];
1499       mean06 +=histo[median+idelta]*(median+idelta);
1500       rms06  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1501     }
1502     if (count09<0.9*count1){
1503       count09+=histo[median-idelta];
1504       mean09 +=histo[median-idelta]*(median-idelta);
1505       rms09  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1506       count09+=histo[median+idelta];
1507       mean09 +=histo[median+idelta]*(median+idelta);
1508       rms09  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1509     }
1510     if (count10<0.95*count1){
1511       count10+=histo[median-idelta];
1512       mean +=histo[median-idelta]*(median-idelta);
1513       rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1514       count10+=histo[median+idelta];
1515       mean +=histo[median+idelta]*(median+idelta);
1516       rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1517     }
1518   }
1519   if (count10) {
1520     mean  /=count10;
1521     rms    = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1522   }
1523   if (count06) {
1524     mean06/=count06;
1525     rms06    = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1526   }
1527   if (count09) {
1528     mean09/=count09;
1529     rms09    = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1530   }
1531   rmsEvent = rms09;
1532   //
1533   pedestalEvent = median;
1534   if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1535   //
1536   UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1537   //
1538   // Dump mean signal info
1539   //
1540     if (AliTPCReconstructor::StreamLevel()>0) {
1541     (*fDebugStreamer)<<"Signal"<<
1542     "TimeStamp="<<fTimeStamp<<
1543     "EventType="<<fEventType<<
1544     "Sector="<<uid[0]<<
1545     "Row="<<uid[1]<<
1546     "Pad="<<uid[2]<<
1547     "Max="<<max<<
1548     "MaxPos="<<maxPos<<
1549     //
1550     "Median="<<median<<
1551     "Mean="<<mean<<
1552     "RMS="<<rms<<
1553     "Mean06="<<mean06<<
1554     "RMS06="<<rms06<<
1555     "Mean09="<<mean09<<
1556     "RMS09="<<rms09<<
1557     "RMSCalib="<<rmsCalib<<
1558     "PedCalib="<<pedestalCalib<<
1559     "\n";
1560     }
1561   //
1562   // fill pedestal histogram
1563   //
1564   //
1565   //
1566   //
1567   Float_t kMin =fRecoParam->GetDumpAmplitudeMin();   // minimal signal to be dumped
1568   Float_t *dsignal = new Float_t[nchannels];
1569   Float_t *dtime   = new Float_t[nchannels];
1570   for (Int_t i=0; i<nchannels; i++){
1571     dtime[i] = i;
1572     dsignal[i] = signal[i];
1573   }
1574
1575   //
1576   // Big signals dumping
1577   //
1578   if (AliTPCReconstructor::StreamLevel()>0) {
1579   if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin())
1580     (*fDebugStreamer)<<"SignalB"<<     // pads with signal
1581       "TimeStamp="<<fTimeStamp<<
1582       "EventType="<<fEventType<<
1583       "Sector="<<uid[0]<<
1584       "Row="<<uid[1]<<
1585       "Pad="<<uid[2]<<
1586       //      "Graph="<<graph<<
1587       "Max="<<max<<
1588       "MaxPos="<<maxPos<<
1589       //
1590       "Median="<<median<<
1591       "Mean="<<mean<<
1592       "RMS="<<rms<<
1593       "Mean06="<<mean06<<
1594       "RMS06="<<rms06<<
1595       "Mean09="<<mean09<<
1596       "RMS09="<<rms09<<
1597       "\n";
1598   }
1599
1600   delete [] dsignal;
1601   delete [] dtime;
1602   if (rms06>fRecoParam->GetMaxNoise()) {
1603     pedestalEvent+=1024.;
1604     return 1024+median; // sign noisy channel in debug mode
1605   }
1606   return median;
1607 }
1608
1609 Int_t AliTPCclustererMI::ReadHLTClusters()
1610 {
1611   //
1612   // read HLT clusters instead of off line custers,
1613   // used in Digits2Clusters
1614   //
1615
1616   TObject* pClusterAccess=NULL;
1617   TClass* pCl=NULL;
1618   ROOT::NewFunc_t pNewFunc=NULL;
1619   do {
1620     pCl=TClass::GetClass("AliHLTTPCClusterAccessHLTOUT");
1621   } while (!pCl && gSystem->Load("libAliHLTTPC.so")==0);
1622   if (!pCl || (pNewFunc=pCl->GetNew())==NULL) {
1623     AliError("can not load class description of AliHLTTPCClusterAccessHLTOUT, aborting ...");
1624     return -1;
1625   }
1626
1627   void* p=(*pNewFunc)(NULL);
1628   if (!p) {
1629     AliError("unable to create instance of AliHLTTPCClusterAccessHLTOUT");
1630     return -2;
1631   }
1632   pClusterAccess=reinterpret_cast<TObject*>(p);
1633   if (!pClusterAccess) {
1634     AliError("instance not of type TObject");
1635     return -3 ;
1636   }
1637
1638   const Int_t kNIS = fParam->GetNInnerSector();
1639   const Int_t kNOS = fParam->GetNOuterSector();
1640   const Int_t kNS = kNIS + kNOS;
1641   fNclusters  = 0;
1642
1643   for(fSector = 0; fSector < kNS; fSector++) {
1644
1645     TString param("sector="); param+=fSector;
1646     pClusterAccess->Clear();
1647     pClusterAccess->Execute("read", param);
1648     if (pClusterAccess->FindObject("clusterarray")==NULL) {
1649       AliError("HLT clusters requested, but not cluster array not present");
1650       return -4;
1651     }
1652
1653     TClonesArray* clusterArray=dynamic_cast<TClonesArray*>(pClusterAccess->FindObject("clusterarray"));
1654     if (!clusterArray) {
1655       AliError("HLT cluster array is not of class type TClonesArray");
1656       return -5;
1657     }
1658
1659     AliDebug(4,Form("Reading %d clusters from HLT for sector %d", clusterArray->GetEntriesFast(), fSector));
1660
1661     Int_t nClusterSector=0;
1662     Int_t nRows=fParam->GetNRow(fSector);
1663
1664     for (fRow = 0; fRow < nRows; fRow++) {
1665       fRowCl->SetID(fParam->GetIndex(fSector, fRow));
1666       if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
1667       fNcluster=0; // reset clusters per row
1668
1669       fRx = fParam->GetPadRowRadii(fSector, fRow);
1670       fPadLength = fParam->GetPadPitchLength(fSector, fRow);
1671       fPadWidth  = fParam->GetPadPitchWidth();
1672       fMaxPad = fParam->GetNPads(fSector,fRow);
1673       fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
1674
1675       fBins = fAllBins[fRow];
1676       fSigBins = fAllSigBins[fRow];
1677       fNSigBins = fAllNSigBins[fRow];
1678
1679       for (Int_t i=0; i<clusterArray->GetEntriesFast(); i++) {
1680         if (!clusterArray->At(i))
1681           continue;
1682
1683         AliTPCclusterMI* cluster=dynamic_cast<AliTPCclusterMI*>(clusterArray->At(i));
1684         if (!cluster) continue;
1685         if (cluster->GetRow()!=fRow) continue;
1686         nClusterSector++;
1687         AddCluster(*cluster, NULL, 0);
1688       }
1689
1690       FillRow();
1691       fRowCl->GetArray()->Clear("c");
1692
1693     } // for (fRow = 0; fRow < nRows; fRow++) {
1694     if (nClusterSector!=clusterArray->GetEntriesFast()) {
1695       AliError(Form("Failed to read %d out of %d HLT clusters",
1696                     clusterArray->GetEntriesFast()-nClusterSector,
1697                     clusterArray->GetEntriesFast()));
1698     }
1699     fNclusters+=nClusterSector;
1700   } // for(fSector = 0; fSector < kNS; fSector++) {
1701
1702   delete pClusterAccess;
1703
1704   Info("Digits2Clusters", "Number of converted HLT clusters : %d", fNclusters);
1705
1706   return 0;
1707 }