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