TPC/AliTPCclustererMI.cxx

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