Some of the coding violations corrected
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.cxx
index f721615..82bb79a 100644 (file)
 //-------------------------------------------------------
 //          Implementation of the TPC clusterer
 //
+//  1. The Input data for reconstruction - Options
+//      1.a Simulated data  - TTree - invoked Digits2Clusters()
+//      1.b Raw data        - Digits2Clusters(AliRawReader* rawReader); 
+//
+//  2. The Output data
+//      2.a TTree with clusters - if  SetOutput(TTree * tree) invoked
+//      2.b TObjArray           - Faster option for HLT
+//      2.c TClonesArray        - Faster option for HLT (smaller memory consumption), activate with fBClonesArray flag
+//
+//  3. Reconstruction setup
+//     see AliTPCRecoParam for list of parameters 
+//     The reconstruction parameterization taken from the 
+//     AliTPCReconstructor::GetRecoParam()
+//     Possible to setup it in reconstruction macro  AliTPCReconstructor::SetRecoParam(...)
+//     
+//
+//
 //   Origin: Marian Ivanov 
 //-------------------------------------------------------
 
-#include "AliTPCclustererMI.h"
-#include "AliTPCclusterMI.h"
-#include <TObjArray.h>
+#include "Riostream.h"
+#include <TF1.h>
 #include <TFile.h>
-#include "AliTPCClustersArray.h"
-#include "AliTPCClustersRow.h"
+#include <TGraph.h>
+#include <TH1F.h>
+#include <TObjArray.h>
+#include <TClonesArray.h>
+#include <TRandom.h>
+#include <TTree.h>
+#include <TTreeStream.h>
+
 #include "AliDigits.h"
+#include "AliLoader.h"
+#include "AliLog.h"
+#include "AliMathBase.h"
+#include "AliRawEventHeaderBase.h"
+#include "AliRawReader.h"
+#include "AliRunLoader.h"
 #include "AliSimDigits.h"
+#include "AliTPCCalPad.h"
+#include "AliTPCCalROC.h"
+#include "AliTPCClustersArray.h"
+#include "AliTPCClustersRow.h"
 #include "AliTPCParam.h"
-#include "Riostream.h"
-#include <TTree.h>
+#include "AliTPCRawStream.h"
+#include "AliTPCRawStreamV3.h"
+#include "AliTPCRecoParam.h"
+#include "AliTPCReconstructor.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCclusterInfo.h"
+#include "AliTPCclusterMI.h"
+#include "AliTPCTransform.h"
+#include "AliTPCclustererMI.h"
 
 ClassImp(AliTPCclustererMI)
 
 
 
-AliTPCclustererMI::AliTPCclustererMI()
+AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam):
+  fBins(0),
+  fSigBins(0),
+  fNSigBins(0),
+  fLoop(0),
+  fMaxBin(0),
+  fMaxTime(1006), // 1000>940 so use 1000, add 3 virtual time bins before and 3 after
+  fMaxPad(0),
+  fSector(-1),
+  fRow(-1),
+  fSign(0),
+  fRx(0),
+  fPadWidth(0),
+  fPadLength(0),
+  fZWidth(0),
+  fPedSubtraction(kFALSE),
+  fEventHeader(0),
+  fTimeStamp(0),
+  fEventType(0),
+  fInput(0),
+  fOutput(0),
+  fOutputArray(0),
+  fOutputClonesArray(0),
+  fRowCl(0),
+  fRowDig(0),
+  fParam(0),
+  fNcluster(0),
+  fNclusters(0),
+  fDebugStreamer(0),
+  fRecoParam(0),
+  fBDumpSignal(kFALSE),
+  fBClonesArray(kFALSE),
+  fAllBins(NULL),
+  fAllSigBins(NULL),
+  fAllNSigBins(NULL)
 {
+  //
+  // COSNTRUCTOR
+  // param     - tpc parameters for given file
+  // recoparam - reconstruction parameters 
+  //
   fInput =0;
-  fOutput=0;
+  fParam = par;
+  if (recoParam) {
+    fRecoParam = recoParam;
+  }else{
+    //set default parameters if not specified
+    fRecoParam = AliTPCReconstructor::GetRecoParam();
+    if (!fRecoParam)  fRecoParam = AliTPCRecoParam::GetLowFluxParam();
+  }
+  if(AliTPCReconstructor::StreamLevel()>0) {
+    fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
+  }
+
+  //  Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
+  fRowCl= new AliTPCClustersRow("AliTPCclusterMI");
+
+  // Non-persistent arrays
+  //
+  //alocate memory for sector - maximal case
+  //
+  AliTPCROC * roc = AliTPCROC::Instance();
+  Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
+  Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
+
+  fAllBins = new Float_t*[nRowsMax];
+  fAllSigBins = new Int_t*[nRowsMax];
+  fAllNSigBins = new Int_t[nRowsMax];
+  for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
+    //
+    Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
+    fAllBins[iRow] = new Float_t[maxBin];
+    memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
+    fAllSigBins[iRow] = new Int_t[maxBin];
+    fAllNSigBins[iRow]=0;
+  }
+}
+//______________________________________________________________
+AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
+              :TObject(param),
+  fBins(0),
+  fSigBins(0),
+  fNSigBins(0),
+  fLoop(0),
+  fMaxBin(0),
+  fMaxTime(0),
+  fMaxPad(0),
+  fSector(-1),
+  fRow(-1),
+  fSign(0),
+  fRx(0),
+  fPadWidth(0),
+  fPadLength(0),
+  fZWidth(0),
+  fPedSubtraction(kFALSE),
+  fEventHeader(0),
+  fTimeStamp(0),
+  fEventType(0),
+  fInput(0),
+  fOutput(0),
+  fOutputArray(0),
+  fOutputClonesArray(0),
+  fRowCl(0),
+  fRowDig(0),
+  fParam(0),
+  fNcluster(0),
+  fNclusters(0),
+  fDebugStreamer(0),
+  fRecoParam(0),
+  fBDumpSignal(kFALSE),
+  fBClonesArray(kFALSE),
+  fAllBins(NULL),
+  fAllSigBins(NULL),
+  fAllNSigBins(NULL)
+{
+  //
+  // dummy
+  //
+  fMaxBin = param.fMaxBin;
+}
+//______________________________________________________________
+AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & param)
+{
+  //
+  // assignment operator - dummy
+  //
+  fMaxBin=param.fMaxBin;
+  return (*this);
+}
+//______________________________________________________________
+AliTPCclustererMI::~AliTPCclustererMI(){
+  //
+  //
+  //
+  if (fDebugStreamer) delete fDebugStreamer;
+  if (fOutputArray){
+    //fOutputArray->Delete();
+    delete fOutputArray;
+  }
+  if (fOutputClonesArray){
+    fOutputClonesArray->Delete();
+    delete fOutputClonesArray;
+  }
+
+  if (fRowCl) {
+    fRowCl->GetArray()->Delete();
+    delete fRowCl;
+  }
+
+  AliTPCROC * roc = AliTPCROC::Instance();
+  Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
+  for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
+    delete [] fAllBins[iRow];
+    delete [] fAllSigBins[iRow];
+  }
+  delete [] fAllBins;
+  delete [] fAllSigBins;
+  delete [] fAllNSigBins;
 }
+
 void AliTPCclustererMI::SetInput(TTree * tree)
 {
   //
@@ -58,19 +253,35 @@ void AliTPCclustererMI::SetInput(TTree * tree)
 void AliTPCclustererMI::SetOutput(TTree * tree) 
 {
   //
+  // Set the output tree
+  // If not set the ObjArray used - Option for HLT 
   //
-  fOutput= tree;  
-  AliTPCClustersRow clrow;
-  AliTPCClustersRow *pclrow=&clrow;  
-  clrow.SetClass("AliTPCclusterMI");
-  clrow.SetArray(1); // to make Clones array
+  if (!tree) return;
+  fOutput= tree;
+  AliTPCClustersRow clrow, *pclrow=&clrow;  
+  pclrow = new AliTPCClustersRow("AliTPCclusterMI");
   fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);    
 }
 
 
+void AliTPCclustererMI::FillRow(){
+  //
+  // fill the output container - 
+  // 2 Options possible
+  //          Tree       
+  //          TObjArray
+  //
+  if (fOutput) fOutput->Fill();
+  if (!fOutput && !fBClonesArray){
+    //
+    if (!fOutputArray) fOutputArray = new TObjArray(fParam->GetNRowsTotal());
+    if (fRowCl && fRowCl->GetArray()->GetEntriesFast()>0) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID());
+  }
+}
+
 Float_t  AliTPCclustererMI::GetSigmaY2(Int_t iz){
   // sigma y2 = in digits  - we don't know the angle
-  Float_t z = iz*fParam->GetZWidth();
+  Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
     (fPadWidth*fPadWidth);
   Float_t sres = 0.25;
@@ -81,9 +292,9 @@ Float_t  AliTPCclustererMI::GetSigmaY2(Int_t iz){
 
 Float_t  AliTPCclustererMI::GetSigmaZ2(Int_t iz){
   //sigma z2 = in digits - angle estimated supposing vertex constraint
-  Float_t z = iz*fZWidth;
+  Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
-  Float_t angular = fPadLength*(fParam->GetZLength()-z)/(fRx*fZWidth);
+  Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
   angular*=angular;
   angular/=12.;
   Float_t sres = fParam->GetZSigma()/fZWidth;
@@ -92,19 +303,23 @@ Float_t  AliTPCclustererMI::GetSigmaZ2(Int_t iz){
   return res;
 }
 
-void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Int_t *bins, UInt_t m,
+void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
 AliTPCclusterMI &c) 
 {
+  //
+  //  k    - Make cluster at position k  
+  //  bins - 2 D array of signals mapped to 1 dimensional array - 
+  //  max  - the number of time bins er one dimension
+  //  c    - refernce to cluster to be filled
+  //
   Int_t i0=k/max;  //central pad
   Int_t j0=k%max;  //central time bin
 
   // set pointers to data
   //Int_t dummy[5] ={0,0,0,0,0};
-  Int_t * matrix[5]; //5x5 matrix with digits  - indexing i = 0 ..4  j = -2..2
-  Int_t * resmatrix[5];
+  Float_t * matrix[5]; //5x5 matrix with digits  - indexing i = 0 ..4  j = -2..2
   for (Int_t di=-2;di<=2;di++){
     matrix[di+2]  =  &bins[k+di*max];
-    resmatrix[di+2]  =  &fResBins[k+di*max];
   }
   //build matrix with virtual charge
   Float_t sigmay2= GetSigmaY2(j0);
@@ -113,7 +328,7 @@ AliTPCclusterMI &c)
   Float_t vmatrix[5][5];
   vmatrix[2][2] = matrix[2][0];
   c.SetType(0);
-  c.SetMax(Short_t(vmatrix[2][2])); // write maximal amplitude
+  c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
   for (Int_t di =-1;di <=1;di++)
     for (Int_t dj =-1;dj <=1;dj++){
       Float_t amp = matrix[di+2][dj];
@@ -181,8 +396,10 @@ AliTPCclusterMI &c)
   
   //
   if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
-  if ( (ry <1.2) && (rz<1.2) ) {
-    //if cluster looks like expected 
+  if ( ((ry <1.2) && (rz<1.2)) || (!fRecoParam->GetDoUnfold())) {
+    //
+    //if cluster looks like expected or Unfolding not switched on
+    //standard COG is used
     //+1.2 deviation from expected sigma accepted
     //    c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
 
@@ -190,26 +407,20 @@ AliTPCclusterMI &c)
     meanj +=j0;
     //set cluster parameters
     c.SetQ(sumw);
-    c.SetY(meani*fPadWidth); 
-    c.SetZ(meanj*fZWidth); 
+    c.SetPad(meani-2.5);
+    c.SetTimeBin(meanj-3);
     c.SetSigmaY2(mi2);
     c.SetSigmaZ2(mj2);
-    AddCluster(c);
-    //remove cluster data from data
-    for (Int_t di=-2;di<=2;di++)
-      for (Int_t dj=-2;dj<=2;dj++){
-       resmatrix[di+2][dj] -= Int_t(vmatrix[di+2][dj+2]);
-       if (resmatrix[di+2][dj]<0) resmatrix[di+2][dj]=0;
-      }
-    resmatrix[2][0] =0;
+    c.SetType(0);
+    AddCluster(c,(Float_t*)vmatrix,k);
     return;     
   }
   //
   //unfolding when neccessary  
   //
   
-  Int_t * matrix2[7]; //7x7 matrix with digits  - indexing i = 0 ..6  j = -3..3
-  Int_t dummy[7]={0,0,0,0,0,0};
+  Float_t * matrix2[7]; //7x7 matrix with digits  - indexing i = 0 ..6  j = -3..3
+  Float_t dummy[7]={0,0,0,0,0,0};
   for (Int_t di=-3;di<=3;di++){
     matrix2[di+3] =  &bins[k+di*max];
     if ((k+di*max)<3)  matrix2[di+3] = &dummy[3];
@@ -225,23 +436,21 @@ AliTPCclusterMI &c)
   meanj +=j0;
   //set cluster parameters
   c.SetQ(sumu);
-  c.SetY(meani*fPadWidth); 
-  c.SetZ(meanj*fZWidth); 
+  c.SetPad(meani-2.5);
+  c.SetTimeBin(meanj-3);
   c.SetSigmaY2(mi2);
   c.SetSigmaZ2(mj2);
   c.SetType(Char_t(overlap)+1);
-  AddCluster(c);
+  AddCluster(c,(Float_t*)vmatrix,k);
 
   //unfolding 2
   meani-=i0;
   meanj-=j0;
-  if (gDebug>4)
-    printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
 }
 
 
 
-void AliTPCclustererMI::UnfoldCluster(Int_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj, 
+void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj, 
                                      Float_t & sumu, Float_t & overlap )
 {
   //
@@ -273,7 +482,7 @@ void AliTPCclustererMI::UnfoldCluster(Int_t * matrix2[7], Float_t recmatrix[5][5
     else{
       Float_t ratio =1;
       if (  ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
-          sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
+           (sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 )){
        Float_t xm2 = sum3i[-dk+3];
        Float_t xm1 = sum3i[+3];
        Float_t x1  = sum3i[2*dk+3];
@@ -364,8 +573,6 @@ void AliTPCclustererMI::UnfoldCluster(Int_t * matrix2[7], Float_t recmatrix[5][5
          }
       }
   }
-  if (gDebug>4) 
-    printf("%f\n", recmatrix[2][2]);
   
 }
 
@@ -388,202 +595,940 @@ Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, F
   return max;
 }
 
-void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){
+void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * /*matrix*/, Int_t /*pos*/){
   //
-  // transform cluster to the global coordinata
-  // add the cluster to the array
   //
-  Float_t meani = c.GetY()/fPadWidth;
-  Float_t meanj = c.GetZ()/fZWidth;
+  // Transform cluster to the rotated global coordinata
+  // Assign labels to the cluster
+  // add the cluster to the array
+  // for more details - See  AliTPCTranform::Transform(x,i,0,1) 
+  Float_t meani = c.GetPad();
+  Float_t meanj = c.GetTimeBin();
 
-  Int_t ki = TMath::Nint(meani-3);
+  Int_t ki = TMath::Nint(meani);
   if (ki<0) ki=0;
   if (ki>=fMaxPad) ki = fMaxPad-1;
-  Int_t kj = TMath::Nint(meanj-3);
+  Int_t kj = TMath::Nint(meanj);
   if (kj<0) kj=0;
   if (kj>=fMaxTime-3) kj=fMaxTime-4;
-  // ki and kj shifted to "real" coordinata
-  c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
-  c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
-  c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
-  
+  // ki and kj shifted as integers coordinata
+  if (fRowDig) {
+    c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
+    c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
+    c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
+  }
   
+  c.SetRow(fRow);
+  c.SetDetector(fSector);
   Float_t s2 = c.GetSigmaY2();
   Float_t w=fParam->GetPadPitchWidth(fSector);
-  
   c.SetSigmaY2(s2*w*w);
   s2 = c.GetSigmaZ2(); 
-  w=fZWidth;
-  c.SetSigmaZ2(s2*w*w);
-  c.SetY((meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
-  c.SetZ(fZWidth*(meanj-3)); 
-  c.SetZ(c.GetZ()   - 3.*fParam->GetZSigma()); // PASA delay 
-  c.SetZ(fSign*(fParam->GetZLength() - c.GetZ()));
+  c.SetSigmaZ2(s2*fZWidth*fZWidth);
+  //
+  //
+  //
   
+  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+  if (!transform) {
+    AliFatal("Tranformations not in calibDB");    
+    return;
+  }
+  transform->SetCurrentRecoParam((AliTPCRecoParam*)fRecoParam);
+  Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()};
+  Int_t i[1]={fSector};
+  transform->Transform(x,i,0,1);
+  c.SetX(x[0]);
+  c.SetY(x[1]);
+  c.SetZ(x[2]);
+  //
+  //
   if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
-    //c.SetSigmaY2(c.GetSigmaY2()*25.);
-    //c.SetSigmaZ2(c.GetSigmaZ2()*4.);
     c.SetType(-(c.GetType()+3));  //edge clusters
   }
   if (fLoop==2) c.SetType(100);
+  if (!AcceptCluster(&c)) return;
 
-  TClonesArray * arr = fRowCl->GetArray();
-  // AliTPCclusterMI * cl = 
-  new ((*arr)[fNcluster]) AliTPCclusterMI(c);
+  // select output 
+  TClonesArray * arr = 0;
+  AliTPCclusterMI * cl = 0;
+
+  if(fBClonesArray==kFALSE) {
+     arr = fRowCl->GetArray();
+     cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
+  } else {
+     cl = new ((*fOutputClonesArray)[fNclusters+fNcluster]) AliTPCclusterMI(c);
+  }
+
+  // if (fRecoParam->DumpSignal() &&matrix ) {
+//     Int_t nbins=0;
+//     Float_t *graph =0;
+//     if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
+//       nbins = fMaxTime;
+//       graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
+//     }
+//     AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
+//     cl->SetInfo(info);
+//   }
+  if (!fRecoParam->DumpSignal()) {
+    cl->SetInfo(0);
+  }
+  
+  if (AliTPCReconstructor::StreamLevel()>1) {
+    Float_t xyz[3];
+    cl->GetGlobalXYZ(xyz);
+     (*fDebugStreamer)<<"Clusters"<<
+       "Cl.="<<cl<<
+       "gx="<<xyz[0]<<
+       "gy="<<xyz[1]<<
+       "gz="<<xyz[2]<<
+       "\n";
+  }
 
   fNcluster++;
 }
 
 
 //_____________________________________________________________________________
-void AliTPCclustererMI::Digits2Clusters(const AliTPCParam *par, Int_t eventn)
+void AliTPCclustererMI::Digits2Clusters()
 {
   //-----------------------------------------------------------------
   // This is a simple cluster finder.
   //-----------------------------------------------------------------
-  TDirectory *savedir=gDirectory; 
 
- if (fInput==0){ 
-    cerr<<"AliTPC::Digits2Clusters(): input tree not initialised !\n";
+  if (!fInput) { 
+    Error("Digits2Clusters", "input tree not initialised");
     return;
   }
-  if (fOutput==0) {
-     cerr<<"AliTPC::Digits2Clusters(): output tree not initialised !\n";
-     return;
+  fRecoParam = AliTPCReconstructor::GetRecoParam();
+  if (!fRecoParam){
+    AliFatal("Can not get the reconstruction parameters");
+  }
+  if(AliTPCReconstructor::StreamLevel()>5) {
+    AliInfo("Parameter Dumps");
+    fParam->Dump();
+    fRecoParam->Dump();
   }
 
+  AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
+  AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
   AliSimDigits digarr, *dummy=&digarr;
   fRowDig = dummy;
   fInput->GetBranch("Segment")->SetAddress(&dummy);
   Stat_t nentries = fInput->GetEntries();
   
-  fMaxTime=par->GetMaxTBin()+6; // add 3 virtual time bins before and 3   after
+  fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3   after
     
-  fParam = par;
-  ((AliTPCParam*)par)->Write(par->GetTitle());
-
   Int_t nclusters  = 0;
-  
+
   for (Int_t n=0; n<nentries; n++) {
     fInput->GetEvent(n);
-    Int_t row;
-    if (!par->AdjustSectorRow(digarr.GetID(),fSector,row)) {
+    if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
       continue;
     }
-        
-    AliTPCClustersRow *clrow= new AliTPCClustersRow();
-    fRowCl = clrow;
-    clrow->SetClass("AliTPCclusterMI");
-    clrow->SetArray(1);
-
-    clrow->SetID(digarr.GetID());
-    fOutput->GetBranch("Segment")->SetAddress(&clrow);
-    fRx=par->GetPadRowRadii(fSector,row);
+    Int_t row = fRow;
+    AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector);  // pad gains per given sector
+    AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector); // noise per given sector
+    //
+
+    fRowCl->SetID(digarr.GetID());
+    if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
+    fRx=fParam->GetPadRowRadii(fSector,row);
     
     
-    const Int_t kNIS=par->GetNInnerSector(), kNOS=par->GetNOuterSector();
+    const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
     fZWidth = fParam->GetZWidth();
     if (fSector < kNIS) {
-      fMaxPad = par->GetNPadsLow(row);
+      fMaxPad = fParam->GetNPadsLow(row);
       fSign = (fSector < kNIS/2) ? 1 : -1;
-      fPadLength = par->GetPadPitchLength(fSector,row);
-      fPadWidth = par->GetPadPitchWidth();
+      fPadLength = fParam->GetPadPitchLength(fSector,row);
+      fPadWidth = fParam->GetPadPitchWidth();
     } else {
-      fMaxPad = par->GetNPadsUp(row);
+      fMaxPad = fParam->GetNPadsUp(row);
       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
-      fPadLength = par->GetPadPitchLength(fSector,row);
-      fPadWidth  = par->GetPadPitchWidth();
+      fPadLength = fParam->GetPadPitchLength(fSector,row);
+      fPadWidth  = fParam->GetPadPitchWidth();
     }
     
     
     fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
-    fBins    =new Int_t[fMaxBin];
-    fResBins =new Int_t[fMaxBin];  //fBins with residuals after 1 finder loop 
-    memset(fBins,0,sizeof(Int_t)*fMaxBin);
+    fBins    =new Float_t[fMaxBin];
+    fSigBins =new Int_t[fMaxBin];
+    fNSigBins = 0;
+    memset(fBins,0,sizeof(Float_t)*fMaxBin);
     
     if (digarr.First()) //MI change
       do {
-       Short_t dig=digarr.CurrentDigit();
-       if (dig<=par->GetZeroSup()) continue;
+       Float_t dig=digarr.CurrentDigit();
+       if (dig<=fParam->GetZeroSup()) continue;
        Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
-       fBins[i*fMaxTime+j]=dig;
+       Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
+       Int_t bin = i*fMaxTime+j;
+       if (gain>0){
+         fBins[bin]=dig/gain;
+       }else{
+         fBins[bin]=0;
+       }
+       fSigBins[fNSigBins++]=bin;
       } while (digarr.Next());
     digarr.ExpandTrackBuffer();
 
-    //add virtual charge at the edge   
-    for (Int_t i=0; i<fMaxTime; i++){
-      Float_t amp1 = fBins[i+3*fMaxTime]; 
-      Float_t amp0 =0;
-      if (amp1>0){
-       Float_t amp2 = fBins[i+4*fMaxTime];
-       if (amp2==0) amp2=0.5;
-       Float_t sigma2 = GetSigmaY2(i);         
-       amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2); 
-       if (gDebug>4) printf("\n%f\n",amp0);
-      }  
-      fBins[i+2*fMaxTime] = Int_t(amp0);    
-      amp0 = 0;
-      amp1 = fBins[(fMaxPad+2)*fMaxTime+i];
-      if (amp1>0){
-       Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime];
-       if (amp2==0) amp2=0.5;
-       Float_t sigma2 = GetSigmaY2(i);         
-       amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);         
-       if (gDebug>4) printf("\n%f\n",amp0);
-      }        
-      fBins[(fMaxPad+3)*fMaxTime+i] = Int_t(amp0);           
+    FindClusters(noiseROC);
+    FillRow();
+    fRowCl->GetArray()->Clear("C");    
+    nclusters+=fNcluster;    
+
+    delete[] fBins;
+    delete[] fSigBins;
+  }  
+  Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
+}
+
+void AliTPCclustererMI::ProcessSectorData(){
+  //
+  // Process the data for the current sector
+  //
+
+  AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
+  AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
+  AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector);  // pedestal per given sector
+  AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector);  // noise per given sector
+  //check the presence of the calibration
+  if (!noiseROC ||!pedestalROC ) {
+    AliError(Form("Missing calibration per sector\t%d\n",fSector));
+    return;
+  }
+  Int_t  nRows=fParam->GetNRow(fSector);
+  Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
+  Int_t zeroSup = fParam->GetZeroSup();
+  //    if (calcPedestal) {
+  if (kFALSE ) {
+    for (Int_t iRow = 0; iRow < nRows; iRow++) {
+      Int_t maxPad = fParam->GetNPads(fSector, iRow);
+      
+      for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
+    //
+    // Temporary fix for data production - !!!! MARIAN
+    // The noise calibration should take mean and RMS - currently the Gaussian fit used
+    // In case of double peak  - the pad should be rejected
+    //
+    // Line mean - if more than given digits over threshold - make a noise calculation
+    // and pedestal substration
+        if (!calcPedestal && fAllBins[iRow][iPad*fMaxTime+0]<50) continue;
+    //
+        if (fAllBins[iRow][iPad*fMaxTime+0] <1 ) continue;  // no data
+        Float_t *p = &fAllBins[iRow][iPad*fMaxTime+3];
+    //Float_t pedestal = TMath::Median(fMaxTime, p);
+        Int_t id[3] = {fSector, iRow, iPad-3};
+    // calib values
+        Double_t rmsCalib=  noiseROC->GetValue(iRow,iPad-3);
+        Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
+        Double_t rmsEvent       = rmsCalib;
+        Double_t pedestalEvent  = pedestalCalib;
+        ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
+        if (rmsEvent<rmsCalib) rmsEvent = rmsCalib;   // take worst scenario
+        if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
+        
+    //
+        for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
+          Int_t bin = iPad*fMaxTime+iTimeBin;
+          fAllBins[iRow][bin] -= pedestalEvent;
+          if (iTimeBin < fRecoParam->GetFirstBin())
+            fAllBins[iRow][bin] = 0;
+          if (iTimeBin > fRecoParam->GetLastBin())
+            fAllBins[iRow][bin] = 0;
+          if (fAllBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
+            fAllBins[iRow][bin] = 0;
+          if (fAllBins[iRow][bin] < 3.0*rmsEvent)   // 3 sigma cut on RMS
+            fAllBins[iRow][bin] = 0;
+          if (fAllBins[iRow][bin]) fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
+        }
+      }
     }
+  }
+  
+  if (AliTPCReconstructor::StreamLevel()>5) {
+    for (Int_t iRow = 0; iRow < nRows; iRow++) {
+      Int_t maxPad = fParam->GetNPads(fSector,iRow);
+      
+      for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
+        for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
+          Int_t bin = iPad*fMaxTime+iTimeBin;
+          Float_t signal = fAllBins[iRow][bin];
+          if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
+            Double_t x[]={iRow,iPad-3,iTimeBin-3};
+            Int_t i[]={fSector};
+            AliTPCTransform trafo;
+            trafo.Transform(x,i,0,1);
+            Double_t gx[3]={x[0],x[1],x[2]};
+            trafo.RotatedGlobal2Global(fSector,gx);
+        //        fAllSigBins[iRow][fAllNSigBins[iRow]++]
+            Int_t rowsigBins = fAllNSigBins[iRow];
+            Int_t first=fAllSigBins[iRow][0];
+            Int_t last= 0;
+        //        if (rowsigBins>0) fAllSigBins[iRow][fAllNSigBins[iRow]-1];
+            
+            if (AliTPCReconstructor::StreamLevel()>5) {
+              (*fDebugStreamer)<<"Digits"<<
+                "sec="<<fSector<<
+                "row="<<iRow<<
+                "pad="<<iPad<<
+                "time="<<iTimeBin<<
+                "sig="<<signal<<
+                "x="<<x[0]<<
+                "y="<<x[1]<<
+                "z="<<x[2]<<
+                "gx="<<gx[0]<<
+                "gy="<<gx[1]<<
+                "gz="<<gx[2]<<
+    //
+                "rowsigBins="<<rowsigBins<<
+                "first="<<first<<
+                "last="<<last<<
+                "\n";
+            }
+          }
+        }
+      }
+    }
+  }
+  
+    // Now loop over rows and find clusters
+  for (fRow = 0; fRow < nRows; fRow++) {
+    fRowCl->SetID(fParam->GetIndex(fSector, fRow));
+    if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
+    
+    fRx = fParam->GetPadRowRadii(fSector, fRow);
+    fPadLength = fParam->GetPadPitchLength(fSector, fRow);
+    fPadWidth  = fParam->GetPadPitchWidth();
+    fMaxPad = fParam->GetNPads(fSector,fRow);
+    fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
+    
+    fBins = fAllBins[fRow];
+    fSigBins = fAllSigBins[fRow];
+    fNSigBins = fAllNSigBins[fRow];
+    
+    FindClusters(noiseROC);
+    
+    FillRow();
+    if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear("C");
+    fNclusters += fNcluster;
+    
+  } // End of loop to find clusters
+}
+
+
+void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
+{
+//-----------------------------------------------------------------
+// This is a cluster finder for the TPC raw data.
+// The method assumes NO ordering of the altro channels.
+// The pedestal subtraction can be switched on and off
+// using an option of the TPC reconstructor
+//-----------------------------------------------------------------
+  fRecoParam = AliTPCReconstructor::GetRecoParam();
+  if (!fRecoParam){
+    AliFatal("Can not get the reconstruction parameters");
+  }
+  if(AliTPCReconstructor::StreamLevel()>5) {
+    AliInfo("Parameter Dumps");
+    fParam->Dump();
+    fRecoParam->Dump();
+  }
+  fRowDig = NULL;
 
-    memcpy(fResBins,fBins, fMaxBin*2);
+  AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
+  AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
+  //
+  AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
+  fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
+  if (fEventHeader){
+    fTimeStamp = fEventHeader->Get("Timestamp");
+    fEventType = fEventHeader->Get("Type");
+    AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+    transform->SetCurrentTimeStamp(fTimeStamp);
+    transform->SetCurrentRun(rawReader->GetRunNumber());
+  }
+  
+  // creaate one TClonesArray for all clusters
+  if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
+  // reset counter
+  fNclusters  = 0;
+  
+  fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
+//   const Int_t kNIS = fParam->GetNInnerSector();
+//   const Int_t kNOS = fParam->GetNOuterSector();
+//   const Int_t kNS = kNIS + kNOS;
+  fZWidth = fParam->GetZWidth();
+  Int_t zeroSup = fParam->GetZeroSup();
+  //
+  // Clean-up
+  //
+  AliTPCROC * roc = AliTPCROC::Instance();
+  Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
+  Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
+  for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
     //
-    fNcluster=0;
-    //first loop - for "gold cluster" 
-    fLoop=1;
-    Int_t *b=&fBins[-1]+2*fMaxTime;
-    Int_t crtime = Int_t((fParam->GetZLength()-1.05*fRx)/fZWidth-5);
-
-    for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
-      b++;
-      if (*b<8) continue;   //threshold form maxima
-      if (i%fMaxTime<crtime) {
-       Int_t delta = -(i%fMaxTime)+crtime;
-       b+=delta;
-       i+=delta;
-       continue; 
+    Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
+    memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
+    fAllNSigBins[iRow]=0;
+  }
+
+  Int_t prevSector=-1;
+  rawReader->Reset();
+  Int_t digCounter=0;
+  //
+  // Loop over DDLs
+  //
+  const Int_t kNIS = fParam->GetNInnerSector();
+  const Int_t kNOS = fParam->GetNOuterSector();
+  const Int_t kNS = kNIS + kNOS;
+  
+  for(fSector = 0; fSector < kNS; fSector++) {
+    
+    Int_t nRows = 0;
+    Int_t nDDLs = 0, indexDDL = 0;
+    if (fSector < kNIS) {
+      nRows = fParam->GetNRowLow();
+      fSign = (fSector < kNIS/2) ? 1 : -1;
+      nDDLs = 2;
+      indexDDL = fSector * 2;
+    }
+    else {
+      nRows = fParam->GetNRowUp();
+      fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
+      nDDLs = 4;
+      indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
+    }
+    
+    // load the raw data for corresponding DDLs
+    rawReader->Reset();
+    rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
+    
+  while (input.NextDDL()){
+    if (input.GetSector() != fSector)
+      AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
+    
+    //Int_t nRows = fParam->GetNRow(fSector);
+    
+    AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
+    // Begin loop over altro data
+    Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
+    Float_t gain =1;
+    
+    //loop over pads
+    while ( input.NextChannel() ) {
+      Int_t iRow = input.GetRow();
+      if (iRow < 0){
+        continue;
+      }
+      if (iRow >= nRows){
+        AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
+                      iRow, 0, nRows -1));
+        continue;
+      }
+      //pad
+      Int_t iPad = input.GetPad();
+      if (iPad < 0 || iPad >= nPadsMax) {
+        AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
+                      iPad, 0, nPadsMax-1));
+        continue;
       }
-     
-      if (!IsMaximum(*b,fMaxTime,b)) continue;
-      AliTPCclusterMI c;
-      Int_t dummy=0;
-      MakeCluster(i, fMaxTime, fBins, dummy,c);
-      //}
+      gain    = gainROC->GetValue(iRow,iPad);
+      iPad+=3;
+
+      //loop over bunches
+      while ( input.NextBunch() ){
+        Int_t  startTbin    = (Int_t)input.GetStartTimeBin();
+        Int_t  bunchlength  = (Int_t)input.GetBunchLength();
+        const UShort_t *sig = input.GetSignals();
+        for (Int_t iTime = 0; iTime<bunchlength; iTime++){
+          Int_t iTimeBin=startTbin-iTime;
+          if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
+            continue;
+            AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
+                          iTimeBin, 0, iTimeBin -1));
+          }
+          iTimeBin+=3;
+          //signal
+          Float_t signal=(Float_t)sig[iTime];
+          if (!calcPedestal && signal <= zeroSup) continue;
+      
+          if (!calcPedestal) {
+            Int_t bin = iPad*fMaxTime+iTimeBin;
+            if (gain>0){
+              fAllBins[iRow][bin] = signal/gain;
+            }else{
+              fAllBins[iRow][bin] =0;
+            }
+            fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
+          }else{
+            fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
+          }
+          fAllBins[iRow][iPad*fMaxTime+0]+=1.;  // pad with signal
+          
+          // Temporary
+          digCounter++;
+        }// end loop signals in bunch
+      }// end loop bunches
+    } // end loop pads
+    //
+    //
+    //
+    //
+    // Now loop over rows and perform pedestal subtraction
+    if (digCounter==0) continue;
+  } // End of loop over sectors
+  //process last sector
+  if ( digCounter>0 ){
+    ProcessSectorData();
+    for (Int_t iRow = 0; iRow < fParam->GetNRow(fSector); iRow++) {
+      Int_t maxPad = fParam->GetNPads(fSector,iRow);
+      Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
+      memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
+      fAllNSigBins[iRow] = 0;
+    }
+    prevSector=fSector;
+    digCounter=0;
+  }
+  }
+  
+  if (rawReader->GetEventId() && fOutput ){
+    Info("Digits2Clusters", "File  %s Event\t%u\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
+  }
+  
+  if(rawReader->GetEventId()) {
+    Info("Digits2Clusters", "Event\t%u\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
+  }
+  
+  if(fBClonesArray) {
+    //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
+  }
+}
+
+
+
+
+
+void AliTPCclustererMI::Digits2ClustersOld
+(AliRawReader* rawReader)
+{
+//-----------------------------------------------------------------
+// This is a cluster finder for the TPC raw data.
+// The method assumes NO ordering of the altro channels.
+// The pedestal subtraction can be switched on and off
+// using an option of the TPC reconstructor
+//-----------------------------------------------------------------
+  fRecoParam = AliTPCReconstructor::GetRecoParam();
+  if (!fRecoParam){
+    AliFatal("Can not get the reconstruction parameters");
+  }
+  if(AliTPCReconstructor::StreamLevel()>5) {
+    AliInfo("Parameter Dumps");
+    fParam->Dump();
+    fRecoParam->Dump();
+  }
+  fRowDig = NULL;
+
+  AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
+  AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
+  //
+  AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
+  fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
+  if (fEventHeader){
+    fTimeStamp = fEventHeader->Get("Timestamp");  
+    fEventType = fEventHeader->Get("Type");  
+  }
+
+  // creaate one TClonesArray for all clusters
+  if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
+  // reset counter
+  fNclusters  = 0;
+  
+  fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
+  const Int_t kNIS = fParam->GetNInnerSector();
+  const Int_t kNOS = fParam->GetNOuterSector();
+  const Int_t kNS = kNIS + kNOS;
+  fZWidth = fParam->GetZWidth();
+  Int_t zeroSup = fParam->GetZeroSup();
+  //
+  // Clean-up
+  //
+
+  AliTPCROC * roc = AliTPCROC::Instance();
+  Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
+  Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
+  for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
+    //
+    Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
+    memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
+    fAllNSigBins[iRow]=0;
+  }
+  //
+  // Loop over sectors
+  //
+  for(fSector = 0; fSector < kNS; fSector++) {
+
+    Int_t nRows = 0;
+    Int_t nDDLs = 0, indexDDL = 0;
+    if (fSector < kNIS) {
+      nRows = fParam->GetNRowLow();
+      fSign = (fSector < kNIS/2) ? 1 : -1;
+      nDDLs = 2;
+      indexDDL = fSector * 2;
     }
-    //memcpy(fBins,fResBins, fMaxBin*2);
-    //second  loop - for rest cluster 
-    /*        
-    fLoop=2;
-    b=&fResBins[-1]+2*fMaxTime;
-    for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
-      b++;
-      if (*b<25) continue;   // bigger threshold for maxima
-      if (!IsMaximum(*b,fMaxTime,b)) continue;
-      AliTPCclusterMI c;
-      Int_t dummy;
-      MakeCluster(i, fMaxTime, fResBins, dummy,c);
-      //}
+    else {
+      nRows = fParam->GetNRowUp();
+      fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
+      nDDLs = 4;
+      indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
     }
-    */
 
-    fOutput->Fill();
-    delete clrow;    
-    nclusters+=fNcluster;    
-    delete[] fBins;      
-    delete[] fResBins;  
-  }  
-  cerr<<"Number of found clusters : "<<nclusters<<"                        \n";
+    // load the raw data for corresponding DDLs
+    rawReader->Reset();
+    rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
+
+    // select only good sector 
+    if (!input.Next()) continue;
+    if(input.GetSector() != fSector) continue;
+
+    AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
+    
+    for (Int_t iRow = 0; iRow < nRows; iRow++) {
+      Int_t maxPad;
+      if (fSector < kNIS)
+        maxPad = fParam->GetNPadsLow(iRow);
+      else
+        maxPad = fParam->GetNPadsUp(iRow);
+      
+      Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
+      memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
+      fAllNSigBins[iRow] = 0;
+    }
+    
+    Int_t digCounter=0;
+    // Begin loop over altro data
+    Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
+    Float_t gain =1;
+    Int_t lastPad=-1;
+
+    input.Reset();
+    while (input.Next()) {
+      if (input.GetSector() != fSector)
+        AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
+
+
+      Int_t iRow = input.GetRow();
+      if (iRow < 0){
+        continue;
+      }
+
+      if (iRow < 0 || iRow >= nRows){
+        AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
+                     iRow, 0, nRows -1));
+        continue;
+      }
+      //pad
+      Int_t iPad = input.GetPad();
+      if (iPad < 0 || iPad >= nPadsMax) {
+        AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
+                     iPad, 0, nPadsMax-1));
+        continue;
+      }
+      if (iPad!=lastPad){
+        gain    = gainROC->GetValue(iRow,iPad);
+        lastPad = iPad;
+      }
+      iPad+=3;
+      //time
+      Int_t iTimeBin = input.GetTime();
+      if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
+        continue;
+        AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
+                     iTimeBin, 0, iTimeBin -1));
+      }
+      iTimeBin+=3;
+
+      //signal
+      Float_t signal = input.GetSignal();
+      if (!calcPedestal && signal <= zeroSup) continue;
+
+      if (!calcPedestal) {
+        Int_t bin = iPad*fMaxTime+iTimeBin;
+        if (gain>0){
+          fAllBins[iRow][bin] = signal/gain;
+        }else{
+          fAllBins[iRow][bin] =0;
+        }
+        fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin;
+      }else{
+        fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
+      }
+      fAllBins[iRow][iPad*fMaxTime+0]+=1.;  // pad with signal
+
+      // Temporary
+      digCounter++;
+    } // End of the loop over altro data
+    //
+    //
+    //
+    //
+    // Now loop over rows and perform pedestal subtraction
+    if (digCounter==0) continue;
+    ProcessSectorData();
+  } // End of loop over sectors
+
+  if (rawReader->GetEventId() && fOutput ){
+    Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
+  }
   
-  fOutput->Write();
-  savedir->cd();
+  if(rawReader->GetEventId()) {
+    Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
+  }
+
+  if(fBClonesArray) {
+    //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
+  }
 }
+
+void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
+{
+  
+  //
+  // add virtual charge at the edge   
+  //
+  Double_t kMaxDumpSize = 500000;
+  if (!fOutput) {
+    fBDumpSignal =kFALSE;
+  }else{
+    if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE;   //dump signal flag
+  }
+   
+  fNcluster=0;
+  fLoop=1;
+  Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
+  Float_t minMaxCutAbs       = fRecoParam->GetMinMaxCutAbs();
+  Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
+  Float_t minUpDownCutAbs    = fRecoParam->GetMinUpDownCutAbs();
+  Float_t minMaxCutSigma       = fRecoParam->GetMinMaxCutSigma();
+  Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
+  Float_t minUpDownCutSigma    = fRecoParam->GetMinUpDownCutSigma();
+  Int_t   useOnePadCluster     = fRecoParam->GetUseOnePadCluster();
+  for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
+    Int_t i = fSigBins[iSig];
+    if (i%fMaxTime<=crtime) continue;
+    Float_t *b = &fBins[i];
+    //absolute custs
+    if (b[0]<minMaxCutAbs) continue;   //threshold for maxima  
+    //
+    if (useOnePadCluster==0){
+      if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue;  // cut on isolated clusters 
+      if (b[-1]+b[1]<=0) continue;               // cut on isolated clusters
+      if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
+    }
+    //
+    if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue;   //threshold for up down  (TRF) 
+    if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue;   //threshold for left right (PRF)    
+    if (!IsMaximum(*b,fMaxTime,b)) continue;
+    //
+    Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
+    if (noise>fRecoParam->GetMaxNoise()) continue;
+    // sigma cuts
+    if (b[0]<minMaxCutSigma*noise) continue;   //threshold form maxima  
+    if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue;   //threshold for up town TRF 
+    if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue;   //threshold for left right (PRF)    
+  
+    AliTPCclusterMI c;   // default cosntruction  without info
+    Int_t dummy=0;
+    MakeCluster(i, fMaxTime, fBins, dummy,c);
+    
+    //}
+  }
+}
+
+Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
+  //
+  // Currently hack to filter digital noise (15.06.2008)
+  // To be parameterized in the AliTPCrecoParam
+  // More inteligent way  to be used in future
+  // Acces to the proper pedestal file needed
+  //
+  if (cl->GetMax()<400) return kTRUE;
+  Double_t ratio = cl->GetQ()/cl->GetMax();
+  if (cl->GetMax()>700){
+    if ((ratio - int(ratio)>0.8)) return kFALSE;
+  }
+  if ((ratio - int(ratio)<0.95)) return kTRUE;
+  return kFALSE;
+}
+
+
+Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
+  //
+  // process signal on given pad - + streaming of additional information in special mode
+  //
+  // id[0] - sector
+  // id[1] - row
+  // id[2] - pad 
+
+  //
+  // ESTIMATE pedestal and the noise
+  // 
+  const Int_t kPedMax = 100;
+  Float_t  max    =  0;
+  Float_t  maxPos =  0;
+  Int_t    median =  -1;
+  Int_t    count0 =  0;
+  Int_t    count1 =  0;
+  Float_t  rmsCalib   = rmsEvent;       // backup initial value ( from calib)
+  Float_t  pedestalCalib = pedestalEvent;// backup initial value ( from calib)
+  Int_t    firstBin = fRecoParam->GetFirstBin();
+  //
+  UShort_t histo[kPedMax];
+  //memset(histo,0,kPedMax*sizeof(UShort_t));
+  for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
+  for (Int_t i=0; i<fMaxTime; i++){
+    if (signal[i]<=0) continue;
+    if (signal[i]>max && i>firstBin) {
+      max = signal[i];
+      maxPos = i;
+    }
+    if (signal[i]>kPedMax-1) continue;
+    histo[int(signal[i]+0.5)]++;
+    count0++;
+  }
+  //
+  for (Int_t i=1; i<kPedMax; i++){
+    if (count1<count0*0.5) median=i;
+    count1+=histo[i];
+  }
+  // truncated mean  
+  //
+  Float_t count10=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
+  Float_t count06=histo[median] ,mean06=histo[median]*median,  rms06=histo[median]*median*median ;
+  Float_t count09=histo[median] ,mean09=histo[median]*median,  rms09=histo[median]*median*median ;
+  //
+  for (Int_t idelta=1; idelta<10; idelta++){
+    if (median-idelta<=0) continue;
+    if (median+idelta>kPedMax) continue;
+    if (count06<0.6*count1){
+      count06+=histo[median-idelta];
+      mean06 +=histo[median-idelta]*(median-idelta);
+      rms06  +=histo[median-idelta]*(median-idelta)*(median-idelta);
+      count06+=histo[median+idelta];
+      mean06 +=histo[median+idelta]*(median+idelta);
+      rms06  +=histo[median+idelta]*(median+idelta)*(median+idelta);
+    }
+    if (count09<0.9*count1){
+      count09+=histo[median-idelta];
+      mean09 +=histo[median-idelta]*(median-idelta);
+      rms09  +=histo[median-idelta]*(median-idelta)*(median-idelta);
+      count09+=histo[median+idelta];
+      mean09 +=histo[median+idelta]*(median+idelta);
+      rms09  +=histo[median+idelta]*(median+idelta)*(median+idelta);
+    }
+    if (count10<0.95*count1){
+      count10+=histo[median-idelta];
+      mean +=histo[median-idelta]*(median-idelta);
+      rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
+      count10+=histo[median+idelta];
+      mean +=histo[median+idelta]*(median+idelta);
+      rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
+    }
+  }
+  if (count10) {
+    mean  /=count10;
+    rms    = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
+  }
+  if (count06) {
+    mean06/=count06;
+    rms06    = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
+  }
+  if (count09) {
+    mean09/=count09;
+    rms09    = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
+  }
+  rmsEvent = rms09;
+  //
+  pedestalEvent = median;
+  if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
+  //
+  UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
+  //
+  // Dump mean signal info
+  //
+    if (AliTPCReconstructor::StreamLevel()>0) {
+    (*fDebugStreamer)<<"Signal"<<
+    "TimeStamp="<<fTimeStamp<<
+    "EventType="<<fEventType<<
+    "Sector="<<uid[0]<<
+    "Row="<<uid[1]<<
+    "Pad="<<uid[2]<<
+    "Max="<<max<<
+    "MaxPos="<<maxPos<<
+    //
+    "Median="<<median<<
+    "Mean="<<mean<<
+    "RMS="<<rms<<      
+    "Mean06="<<mean06<<
+    "RMS06="<<rms06<<
+    "Mean09="<<mean09<<
+    "RMS09="<<rms09<<
+    "RMSCalib="<<rmsCalib<<
+    "PedCalib="<<pedestalCalib<<
+    "\n";
+    }
+  //
+  // fill pedestal histogram
+  //
+  //
+  //
+  //
+  Float_t kMin =fRecoParam->GetDumpAmplitudeMin();   // minimal signal to be dumped
+  Float_t *dsignal = new Float_t[nchannels];
+  Float_t *dtime   = new Float_t[nchannels];
+  for (Int_t i=0; i<nchannels; i++){
+    dtime[i] = i;
+    dsignal[i] = signal[i];
+  }
+
+  //
+  // Big signals dumping
+  //    
+  if (AliTPCReconstructor::StreamLevel()>0) {
+  if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin()) 
+    (*fDebugStreamer)<<"SignalB"<<     // pads with signal
+      "TimeStamp="<<fTimeStamp<<
+      "EventType="<<fEventType<<
+      "Sector="<<uid[0]<<
+      "Row="<<uid[1]<<
+      "Pad="<<uid[2]<<
+      //      "Graph="<<graph<<
+      "Max="<<max<<
+      "MaxPos="<<maxPos<<      
+      //
+      "Median="<<median<<
+      "Mean="<<mean<<
+      "RMS="<<rms<<      
+      "Mean06="<<mean06<<
+      "RMS06="<<rms06<<
+      "Mean09="<<mean09<<
+      "RMS09="<<rms09<<
+      "\n";
+  }
+
+  delete [] dsignal;
+  delete [] dtime;
+  if (rms06>fRecoParam->GetMaxNoise()) {
+    pedestalEvent+=1024.;
+    return 1024+median; // sign noisy channel in debug mode
+  }
+  return median;
+}
+
+
+
+
+