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