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