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