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