Some of the coding violations corrected
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.cxx
index 0371a87..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 <TGraph.h>
 #include <TH1F.h>
 #include <TObjArray.h>
+#include <TClonesArray.h>
 #include <TRandom.h>
 #include <TTree.h>
 #include <TTreeStream.h>
-#include <TVirtualFFT.h>
 
 #include "AliDigits.h"
 #include "AliLoader.h"
@@ -46,6 +63,7 @@
 #include "AliTPCClustersRow.h"
 #include "AliTPCParam.h"
 #include "AliTPCRawStream.h"
+#include "AliTPCRawStreamV3.h"
 #include "AliTPCRecoParam.h"
 #include "AliTPCReconstructor.h"
 #include "AliTPCcalibDB.h"
@@ -64,7 +82,7 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar
   fNSigBins(0),
   fLoop(0),
   fMaxBin(0),
-  fMaxTime(0),
+  fMaxTime(1006), // 1000>940 so use 1000, add 3 virtual time bins before and 3 after
   fMaxPad(0),
   fSector(-1),
   fRow(-1),
@@ -74,30 +92,32 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar
   fPadLength(0),
   fZWidth(0),
   fPedSubtraction(kFALSE),
-  fIsOldRCUFormat(kFALSE),
   fEventHeader(0),
   fTimeStamp(0),
   fEventType(0),
   fInput(0),
   fOutput(0),
+  fOutputArray(0),
+  fOutputClonesArray(0),
   fRowCl(0),
   fRowDig(0),
   fParam(0),
   fNcluster(0),
-  fAmplitudeHisto(0),
+  fNclusters(0),
   fDebugStreamer(0),
   fRecoParam(0),
   fBDumpSignal(kFALSE),
-  fFFTr2c(0)
+  fBClonesArray(kFALSE),
+  fAllBins(NULL),
+  fAllSigBins(NULL),
+  fAllNSigBins(NULL)
 {
   //
   // COSNTRUCTOR
   // param     - tpc parameters for given file
   // recoparam - reconstruction parameters 
   //
-  fIsOldRCUFormat = kFALSE;
   fInput =0;
-  fOutput=0;
   fParam = par;
   if (recoParam) {
     fRecoParam = recoParam;
@@ -106,10 +126,33 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar
     fRecoParam = AliTPCReconstructor::GetRecoParam();
     if (!fRecoParam)  fRecoParam = AliTPCRecoParam::GetLowFluxParam();
   }
-  fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
-  fAmplitudeHisto = 0;
-  Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
-  fFFTr2c = TVirtualFFT::FFT(1, &nPoints, "R2C  K");
+  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)
@@ -129,21 +172,25 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
   fPadLength(0),
   fZWidth(0),
   fPedSubtraction(kFALSE),
-  fIsOldRCUFormat(kFALSE),
   fEventHeader(0),
   fTimeStamp(0),
   fEventType(0),
   fInput(0),
   fOutput(0),
+  fOutputArray(0),
+  fOutputClonesArray(0),
   fRowCl(0),
   fRowDig(0),
   fParam(0),
   fNcluster(0),
-  fAmplitudeHisto(0),
+  fNclusters(0),
   fDebugStreamer(0),
   fRecoParam(0),
   fBDumpSignal(kFALSE),
-  fFFTr2c(0)
+  fBClonesArray(kFALSE),
+  fAllBins(NULL),
+  fAllSigBins(NULL),
+  fAllNSigBins(NULL)
 {
   //
   // dummy
@@ -161,9 +208,33 @@ AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & para
 }
 //______________________________________________________________
 AliTPCclustererMI::~AliTPCclustererMI(){
-  DumpHistos();
-  if (fAmplitudeHisto) delete fAmplitudeHisto;
+  //
+  //
+  //
   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)
@@ -182,16 +253,32 @@ 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()+fParam->GetNTBinsL1()*fParam->GetZWidth();
@@ -309,7 +396,7 @@ AliTPCclusterMI &c)
   
   //
   if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
-  if ( (ry <1.2) && (rz<1.2) || (!fRecoParam->GetDoUnfold())) {
+  if ( ((ry <1.2) && (rz<1.2)) || (!fRecoParam->GetDoUnfold())) {
     //
     //if cluster looks like expected or Unfolding not switched on
     //standard COG is used
@@ -321,7 +408,7 @@ AliTPCclusterMI &c)
     //set cluster parameters
     c.SetQ(sumw);
     c.SetPad(meani-2.5);
-    c.SetTimeBin(meanj-2.5);
+    c.SetTimeBin(meanj-3);
     c.SetSigmaY2(mi2);
     c.SetSigmaZ2(mj2);
     c.SetType(0);
@@ -359,8 +446,6 @@ AliTPCclusterMI &c)
   //unfolding 2
   meani-=i0;
   meanj-=j0;
-  if (gDebug>4)
-    printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
 }
 
 
@@ -397,7 +482,7 @@ void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[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];
@@ -488,8 +573,6 @@ void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5]
          }
       }
   }
-  if (gDebug>4) 
-    printf("%f\n", recmatrix[2][2]);
   
 }
 
@@ -512,7 +595,7 @@ Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, F
   return max;
 }
 
-void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){
+void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * /*matrix*/, Int_t /*pos*/){
   //
   //
   // Transform cluster to the rotated global coordinata
@@ -545,10 +628,13 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t p
   //
   //
   //
+  
   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
   if (!transform) {
-    AliFatal("Tranformations not in calibDB");
+    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);
@@ -557,32 +643,47 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t p
   c.SetZ(x[2]);
   //
   //
-  if (!fRecoParam->GetBYMirror()){
-    if (fSector%36>17){
-      c.SetY(-c.GetY());
-    }
-  }
-
   if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
     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);
-  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);
+  // 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++;
 }
@@ -599,10 +700,14 @@ void AliTPCclustererMI::Digits2Clusters()
     Error("Digits2Clusters", "input tree not initialised");
     return;
   }
-  if (!fOutput) {
-    Error("Digits2Clusters", "output tree not initialised");
-    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();
@@ -626,13 +731,9 @@ void AliTPCclustererMI::Digits2Clusters()
     AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector);  // pad gains per given sector
     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector); // noise per given sector
     //
-    AliTPCClustersRow *clrow= new AliTPCClustersRow();
-    fRowCl = clrow;
-    clrow->SetClass("AliTPCclusterMI");
-    clrow->SetArray(1);
 
-    clrow->SetID(digarr.GetID());
-    fOutput->GetBranch("Segment")->SetAddress(&clrow);
+    fRowCl->SetID(digarr.GetID());
+    if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
     fRx=fParam->GetPadRowRadii(fSector,row);
     
     
@@ -664,23 +765,161 @@ void AliTPCclustererMI::Digits2Clusters()
        Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
        Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
        Int_t bin = i*fMaxTime+j;
-       fBins[bin]=dig/gain;
+       if (gain>0){
+         fBins[bin]=dig/gain;
+       }else{
+         fBins[bin]=0;
+       }
        fSigBins[fNSigBins++]=bin;
       } while (digarr.Next());
     digarr.ExpandTrackBuffer();
 
     FindClusters(noiseROC);
-
-    fOutput->Fill();
-    delete clrow;    
+    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)
 {
 //-----------------------------------------------------------------
@@ -689,17 +928,211 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
 // 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;
 
-  if (!fOutput) {
-    Error("Digits2Clusters", "output tree not initialised");
-    return;
+  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++) {
+    //
+    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;
+      }
+      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;
-  AliTPCROC * roc = AliTPCROC::Instance();
+
   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
-  AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
-  AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
   //
   AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
@@ -709,8 +1142,10 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
     fEventType = fEventHeader->Get("Type");  
   }
 
-
-  Int_t nclusters  = 0;
+  // 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();
@@ -719,33 +1154,23 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
   fZWidth = fParam->GetZWidth();
   Int_t zeroSup = fParam->GetZeroSup();
   //
-  //alocate memory for sector - maximal case
+  // Clean-up
   //
-  Float_t** allBins = NULL;
-  Int_t** allSigBins = NULL;
-  Int_t*  allNSigBins = NULL;
+
+  AliTPCROC * roc = AliTPCROC::Instance();
   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
-  allBins = new Float_t*[nRowsMax];
-  allSigBins = new Int_t*[nRowsMax];
-  allNSigBins = 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
-    allBins[iRow] = new Float_t[maxBin];
-    memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
-    allSigBins[iRow] = new Int_t[maxBin];
-    allNSigBins[iRow]=0;
+    memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
+    fAllNSigBins[iRow]=0;
   }
   //
   // Loop over sectors
   //
   for(fSector = 0; fSector < kNS; fSector++) {
 
-    AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
-    AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector);  // pedestal per given sector
-    AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector);  // noise per given sector
     Int_t nRows = 0;
     Int_t nDDLs = 0, indexDDL = 0;
     if (fSector < kNIS) {
@@ -761,157 +1186,111 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
     }
 
+    // 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);
+        maxPad = fParam->GetNPadsLow(iRow);
       else
-       maxPad = fParam->GetNPadsUp(iRow);
+        maxPad = fParam->GetNPadsUp(iRow);
       
       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
-      memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
-      allNSigBins[iRow] = 0;
+      memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
+      fAllNSigBins[iRow] = 0;
     }
     
-    // Loas the raw data for corresponding DDLs
-    rawReader->Reset();
-    input.SetOldRCUFormat(fIsOldRCUFormat);
-    rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
     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()));
+        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) !",
+        AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
                      iRow, 0, nRows -1));
-       continue;
+        continue;
       }
       //pad
       Int_t iPad = input.GetPad();
       if (iPad < 0 || iPad >= nPadsMax) {
-       AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
+        AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
                      iPad, 0, nPadsMax-1));
-       continue;
+        continue;
       }
       if (iPad!=lastPad){
-       gain    = gainROC->GetValue(iRow,iPad);
-       lastPad = iPad;
+        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) !",
+        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 && signal <= zeroSup) continue;
+
       if (!calcPedestal) {
-       Int_t bin = iPad*fMaxTime+iTimeBin;
-       allBins[iRow][bin] = signal/gain;
-       allSigBins[iRow][allNSigBins[iRow]++] = bin;
+        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{
-       allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
+        fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
       }
-      allBins[iRow][iPad*fMaxTime+0]=1.;  // pad with 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;
-    //    if (fPedSubtraction) {
-    if (calcPedestal) {
-      for (Int_t iRow = 0; iRow < nRows; iRow++) {
-       Int_t maxPad;
-       if (fSector < kNIS)
-         maxPad = fParam->GetNPadsLow(iRow);
-       else
-         maxPad = fParam->GetNPadsUp(iRow);
-
-       for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
-         if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue;  // no data
-         Float_t *p = &allBins[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;
-           allBins[iRow][bin] -= pedestalEvent;
-           if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())  
-             allBins[iRow][bin] = 0;
-           if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())  
-             allBins[iRow][bin] = 0;
-           if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
-             allBins[iRow][bin] = 0;
-           if (allBins[iRow][bin] < 3.0*rmsEvent)   // 3 sigma cut on RMS
-             allBins[iRow][bin] = 0;
-           if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
-         }
-       }
-      }
-    }
-    // Now loop over rows and find clusters
-    for (fRow = 0; fRow < nRows; fRow++) {
-      fRowCl = new AliTPCClustersRow;
-      fRowCl->SetClass("AliTPCclusterMI");
-      fRowCl->SetArray(1);
-      fRowCl->SetID(fParam->GetIndex(fSector, fRow));
-      fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
-
-      fRx = fParam->GetPadRowRadii(fSector, fRow);
-      fPadLength = fParam->GetPadPitchLength(fSector, fRow);
-      fPadWidth  = fParam->GetPadPitchWidth();
-      if (fSector < kNIS)
-       fMaxPad = fParam->GetNPadsLow(fRow);
-      else
-       fMaxPad = fParam->GetNPadsUp(fRow);
-      fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
-
-      fBins = allBins[fRow];
-      fSigBins = allSigBins[fRow];
-      fNSigBins = allNSigBins[fRow];
-
-      FindClusters(noiseROC);
-
-      fOutput->Fill();
-      delete fRowCl;    
-      nclusters += fNcluster;    
-    } // End of loop to find clusters
+    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);
+  }
   
-  for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
-    delete [] allBins[iRow];
-    delete [] allSigBins[iRow];
-  }  
-  delete [] allBins;
-  delete [] allSigBins;
-  delete [] allNSigBins;
-  
-  Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
+  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)
@@ -921,8 +1300,12 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
   // add virtual charge at the edge   
   //
   Double_t kMaxDumpSize = 500000;
-  if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE;   //dump signal flag
-  //
+  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);
@@ -932,6 +1315,7 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
   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;
@@ -939,21 +1323,24 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
     //absolute custs
     if (b[0]<minMaxCutAbs) continue;   //threshold for maxima  
     //
-    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 (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(kFALSE);   // default cosntruction  without info
+    AliTPCclusterMI c;   // default cosntruction  without info
     Int_t dummy=0;
     MakeCluster(i, fMaxTime, fBins, dummy,c);
     
@@ -961,6 +1348,22 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
   }
 }
 
+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){
   //
@@ -974,7 +1377,6 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
   // ESTIMATE pedestal and the noise
   // 
   const Int_t kPedMax = 100;
-  Double_t kMaxDebugSize = 5000000.;
   Float_t  max    =  0;
   Float_t  maxPos =  0;
   Int_t    median =  -1;
@@ -982,10 +1384,11 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
   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 = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
+  Int_t    firstBin = fRecoParam->GetFirstBin();
   //
   UShort_t histo[kPedMax];
-  memset(histo,0,kPedMax*sizeof(UShort_t));
+  //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) {
@@ -1056,7 +1459,8 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
   //
   // Dump mean signal info
   //
-  (*fDebugStreamer)<<"Signal"<<
+    if (AliTPCReconstructor::StreamLevel()>0) {
+    (*fDebugStreamer)<<"Signal"<<
     "TimeStamp="<<fTimeStamp<<
     "EventType="<<fEventType<<
     "Sector="<<uid[0]<<
@@ -1075,39 +1479,10 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
     "RMSCalib="<<rmsCalib<<
     "PedCalib="<<pedestalCalib<<
     "\n";
+    }
   //
   // fill pedestal histogram
   //
-  AliTPCROC * roc = AliTPCROC::Instance();
-  if (!fAmplitudeHisto){
-    fAmplitudeHisto = new TObjArray(72);
-  }  
-  //
-  if (uid[0]<roc->GetNSectors() 
-      && uid[1]< roc->GetNRows(uid[0])  && 
-      uid[2] <roc->GetNPads(uid[0], uid[1])){
-    TObjArray  * sectorArray = (TObjArray*)fAmplitudeHisto->UncheckedAt(uid[0]);
-    if (!sectorArray){
-      Int_t npads =roc->GetNChannels(uid[0]);
-      sectorArray = new TObjArray(npads);
-      fAmplitudeHisto->AddAt(sectorArray, uid[0]);
-    }
-    Int_t position =  uid[2]+roc->GetRowIndexes(uid[0])[uid[1]];
-    // TH1F * histo = (TH1F*)sectorArray->UncheckedAt(position);
-//     if (!histo){
-//       char hname[100];
-//       sprintf(hname,"Amp_%d_%d_%d",uid[0],uid[1],uid[2]);
-//       TFile * backup = gFile;
-//       fDebugStreamer->GetFile()->cd();
-//       histo = new TH1F(hname, hname, 100, 5,100);
-//       //histo->SetDirectory(0);     // histogram not connected to directory -(File)
-//       sectorArray->AddAt(histo, position);
-//       if (backup) backup->cd();
-//     }
-//     for (Int_t i=0; i<nchannels; i++){
-//       histo->Fill(signal[i]);
-//     }
-  }
   //
   //
   //
@@ -1118,265 +1493,31 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
     dtime[i] = i;
     dsignal[i] = signal[i];
   }
-  //
-  // Digital noise
-  //
- //  if (max-median>30.*TMath::Max(1.,Double_t(rms06)) &&  (((*fDebugStreamer)<<"SignalDN").GetSize()<kMaxDebugSize)){    
-//     //
-//     //
-//     TGraph * graph =new TGraph(nchannels, dtime, dsignal);
-//     //
-//     //
-//     // jumps left - right
-//     Int_t    njumps0=0;
-//     Double_t deltaT0[2000];
-//     Double_t deltaA0[2000];
-//     Int_t    lastJump0 = fRecoParam->GetFirstBin();
-//     Int_t    njumps1=0;
-//     Double_t deltaT1[2000];
-//     Double_t deltaA1[2000];
-//     Int_t    lastJump1 = fRecoParam->GetFirstBin();
-//     Int_t    njumps2=0;
-//     Double_t deltaT2[2000];
-//     Double_t deltaA2[2000];
-//     Int_t    lastJump2 = fRecoParam->GetFirstBin();
-
-//     for (Int_t itime=fRecoParam->GetFirstBin()+1; itime<fRecoParam->GetLastBin()-1; itime++){
-//       if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06))  && 
-//       TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06))  &&
-//       (dsignal[itime-1]-median<5.*rms06) &&
-//       (dsignal[itime+1]-median<5.*rms06)      
-//       ){
-//     deltaA0[njumps0] = dsignal[itime]-dsignal[itime-1];
-//     deltaT0[njumps0] = itime-lastJump0;
-//     lastJump0 = itime;
-//     njumps0++;
-//       }
-//       if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06)) &&
-//       (dsignal[itime-1]-median<5.*rms06) 
-//       ) {
-//     deltaA1[njumps1] = dsignal[itime]-dsignal[itime-1];
-//     deltaT1[njumps1] = itime-lastJump1;
-//     lastJump1 = itime;
-//     njumps1++;
-//       }
-//       if (TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06)) &&
-//       (dsignal[itime+1]-median<5.*rms06) 
-//       ) {
-//     deltaA2[njumps2] = dsignal[itime]-dsignal[itime+1];
-//     deltaT2[njumps2] = itime-lastJump2;
-//     lastJump2 = itime;
-//     njumps2++;
-//       }
-//     }
-//     //
-//     if (njumps0>0 || njumps1>0 || njumps2>0){
-//       TGraph *graphDN0 = new TGraph(njumps0, deltaT0, deltaA0);
-//       TGraph *graphDN1 = new TGraph(njumps1, deltaT1, deltaA1);
-//       TGraph *graphDN2 = new TGraph(njumps2, deltaT2, deltaA2);
-//       (*fDebugStreamer)<<"SignalDN"<<    //digital - noise pads - or random sample of pads
-//     "TimeStamp="<<fTimeStamp<<
-//     "EventType="<<fEventType<<
-//     "Sector="<<uid[0]<<
-//     "Row="<<uid[1]<<
-//     "Pad="<<uid[2]<<
-//     "Graph="<<graph<<
-//     "Max="<<max<<
-//     "MaxPos="<<maxPos<<
-//     "Graph.="<<graph<<  
-//     "P0GraphDN0.="<<graphDN0<<
-//     "P1GraphDN1.="<<graphDN1<<
-//     "P2GraphDN2.="<<graphDN2<<
-//     //
-//     "Median="<<median<<
-//     "Mean="<<mean<<
-//     "RMS="<<rms<<      
-//     "Mean06="<<mean06<<
-//     "RMS06="<<rms06<<
-//     "Mean09="<<mean09<<
-//     "RMS09="<<rms09<<
-//     "\n";
-//       delete graphDN0;
-//       delete graphDN1;
-//       delete graphDN2;
-//     }
-//     delete graph;
-//   }
-
-  //
-  // NOISE STUDY  Fourier transform
-  //
-  TGraph * graph;
-  Bool_t random = (gRandom->Rndm()<0.0003);
-  if (((*fDebugStreamer)<<"SignalN").GetSize()<kMaxDebugSize)
-    if (max-median>kMin || rms06>1.*fParam->GetZeroSup() || random){
-    graph =new TGraph(nchannels, dtime, dsignal);
-    if (rms06>1.*fParam->GetZeroSup() || random){
-      //Double_t *input, Double_t threshold, Bool_t locMax, Double_t *freq, Double_t *re, Double_t *im, Double_t *mag, Double_t *phi);
-      Float_t * input = &(dsignal[fRecoParam->GetFirstBin()]);
-      Float_t freq[2000], re[2000], im[2000], mag[2000], phi[2000];
-      Int_t npoints = TransformFFT(input, -1,kFALSE, freq, re, im, mag, phi);
-      TGraph *graphMag0 = new TGraph(npoints, freq, mag);
-      TGraph *graphPhi0 = new TGraph(npoints, freq, phi);
-      npoints = TransformFFT(input, 0.5,kTRUE, freq, re, im, mag, phi);
-      TGraph *graphMag1 = new TGraph(npoints, freq, mag);
-      TGraph *graphPhi1 = new TGraph(npoints, freq, phi);
-      
-      (*fDebugStreamer)<<"SignalN"<<    //noise pads - or random sample of pads
-       "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<<
-       // FFT part
-       "Mag0.="<<graphMag0<<
-       "Mag1.="<<graphMag1<<
-       "Phi0.="<<graphPhi0<<
-       "Phi1.="<<graphPhi1<<
-       "\n";
-      delete graphMag0;
-      delete graphMag1;
-      delete graphPhi0;
-      delete graphPhi1;
-    }
-    //
-    // Big signals dumping
-    //
-    
-    if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->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 graph;
-  }
-  
-  //
-  //
-  //  Central Electrode signal analysis  
-  //
-  Float_t ceQmax  =0, ceQsum=0, ceTime=0;
-  Float_t cemean  = mean06, cerms=rms06 ;
-  Int_t    cemaxpos= 0;
-  Float_t ceThreshold=5.*cerms;
-  Float_t ceSumThreshold=8.*cerms;
-  const Int_t    kCemin=5;  // range for the analysis of the ce signal +- channels from the peak
-  const Int_t    kCemax=5;
-  for (Int_t i=nchannels-2; i>nchannels/2; i--){
-    if ( (dsignal[i]-mean06)>ceThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] ){
-      cemaxpos=i;
-      break;
-    }
-  }
-  if (cemaxpos!=0){
-    ceQmax = 0;
-    Int_t cemaxpos2=0;
-    for (Int_t i=cemaxpos-20; i<cemaxpos+5; i++){
-      if (i<0 || i>nchannels-1) continue;
-      Double_t val=dsignal[i]- cemean;
-      if (val>ceQmax){
-       cemaxpos2=i;
-       ceQmax = val;
-      }
-    }
-    cemaxpos = cemaxpos2;
-    for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; i++){
-      if (i>0 && i<nchannels&&dsignal[i]- cemean>0){
-       Double_t val=dsignal[i]- cemean;
-       ceTime+=val*dtime[i];
-       ceQsum+=val;
-       if (val>ceQmax) ceQmax=val;
-      }
-    }
-    if (ceQmax&&ceQsum>ceSumThreshold) {
-      ceTime/=ceQsum;
-      (*fDebugStreamer)<<"Signalce"<<
-       "TimeStamp="<<fTimeStamp<<
-       "EventType="<<fEventType<<
-       "Sector="<<uid[0]<<
-       "Row="<<uid[1]<<
-       "Pad="<<uid[2]<<
-       "Max="<<ceQmax<<
-       "Qsum="<<ceQsum<<
-       "Time="<<ceTime<<
-       "RMS06="<<rms06<<
-       //
-       "\n";
-    }
-  }
-  // end of ce signal analysis
-  //
 
   //
-  //  Gating grid signal analysis  
-  //
-  Double_t ggQmax  =0, ggQsum=0, ggTime=0;
-  Double_t ggmean  = mean06, ggrms=rms06 ;
-  Int_t    ggmaxpos= 0;
-  Double_t ggThreshold=5.*ggrms;
-  Double_t ggSumThreshold=8.*ggrms;
-
-  for (Int_t i=1; i<nchannels/4; i++){
-    if ( (dsignal[i]-mean06)>ggThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] &&
-        (dsignal[i]+dsignal[i+1]+dsignal[i-1]-3*mean06)>ggSumThreshold){
-      ggmaxpos=i;
-      if (dsignal[i-1]>dsignal[i+1]) ggmaxpos=i-1;
-      break;
-    }
-  }
-  if (ggmaxpos!=0){
-      for (Int_t i=ggmaxpos-1; i<ggmaxpos+3; i++){       
-         if (i>0 && i<nchannels && dsignal[i]-ggmean>0){
-             Double_t val=dsignal[i]- ggmean;
-             ggTime+=val*dtime[i];
-             ggQsum+=val;
-             if (val>ggQmax) ggQmax=val;
-         }
-      }
-      if (ggQmax&&ggQsum>ggSumThreshold) {
-         ggTime/=ggQsum;
-         (*fDebugStreamer)<<"Signalgg"<<
-           "TimeStamp="<<fTimeStamp<<
-           "EventType="<<fEventType<<
-             "Sector="<<uid[0]<<
-             "Row="<<uid[1]<<
-             "Pad="<<uid[2]<<
-             "Max="<<ggQmax<<
-             "Qsum="<<ggQsum<<
-             "Time="<<ggTime<<
-             "RMS06="<<rms06<<
-             //
-             "\n";
-      }
+  // 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";
   }
-  // end of gg signal analysis
-      
 
   delete [] dsignal;
   delete [] dtime;
@@ -1389,103 +1530,5 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
 
 
 
-void AliTPCclustererMI::DumpHistos(){
-  //
-  // Dump histogram information
-  //
-  if (!fAmplitudeHisto) return;
-  AliTPCROC * roc = AliTPCROC::Instance();
-  for (UInt_t isector=0; isector<AliTPCROC::Instance()->GetNSectors(); isector++){
-    TObjArray * array = (TObjArray*)fAmplitudeHisto->UncheckedAt(isector);
-    if (!array) continue;
-    for (UInt_t ipad = 0; ipad <(UInt_t)array->GetEntriesFast(); ipad++){
-      TH1F * histo = (TH1F*) array->UncheckedAt(ipad);
-      if (!histo) continue;
-      if (histo->GetEntries()<100) continue;
-      histo->Fit("gaus","q");
-      Float_t mean =  histo->GetMean();
-      Float_t rms  =  histo->GetRMS();
-      Float_t gmean = histo->GetFunction("gaus")->GetParameter(1);
-      Float_t gsigma = histo->GetFunction("gaus")->GetParameter(2);
-      Float_t gmeanErr = histo->GetFunction("gaus")->GetParError(1);
-      Float_t gsigmaErr = histo->GetFunction("gaus")->GetParError(2);
-      Float_t max = histo->GetFunction("gaus")->GetParameter(0);
-
-      // get pad number
-      UInt_t row=0, pad =0;
-      const UInt_t *indexes =roc->GetRowIndexes(isector);
-      for (UInt_t irow=0; irow<roc->GetNRows(isector); irow++){
-       if (indexes[irow]<=ipad){
-         row = irow;
-         pad = ipad-indexes[irow];
-       }
-      }      
-      Int_t rpad = pad - (AliTPCROC::Instance()->GetNPads(isector,row))/2;
-      //
-      (*fDebugStreamer)<<"Fit"<<
-       "TimeStamp="<<fTimeStamp<<
-       "EventType="<<fEventType<<
-       "Sector="<<isector<<
-       "Row="<<row<<
-       "Pad="<<pad<<
-       "RPad="<<rpad<<
-       "Max="<<max<<
-       "Mean="<<mean<<
-       "RMS="<<rms<<      
-       "GMean="<<gmean<<
-       "GSigma="<<gsigma<<
-       "GMeanErr="<<gmeanErr<<
-       "GSigmaErr="<<gsigmaErr<<
-       "\n";
-      if (array->UncheckedAt(ipad)) fDebugStreamer->StoreObject(array->UncheckedAt(ipad));
-    }
-  }
-}
 
 
-
-Int_t  AliTPCclustererMI::TransformFFT(Float_t *input, Float_t threshold, Bool_t locMax, Float_t *freq, Float_t *re, Float_t *im, Float_t *mag, Float_t *phi)
-{
-  //
-  // calculate fourrie transform 
-  // return only frequncies with mag over threshold
-  // if locMax is spectified only freque with local maxima over theshold is returned 
-
-  if (! fFFTr2c) return kFALSE;
-  if (!freq) return kFALSE;
-
-  Int_t current=0;
-  Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
-  Double_t *in = new Double_t[nPoints];
-  Double_t *rfft = new Double_t[nPoints];
-  Double_t *ifft = new Double_t[nPoints];
-  for (Int_t i=0; i<nPoints; i++){in[i]=input[i];}
-  fFFTr2c->SetPoints(in);
-  fFFTr2c->Transform();
-  fFFTr2c->GetPointsComplex(rfft, ifft);
-  for (Int_t i=3; i<nPoints/2-3; i++){
-    Float_t lmag =  TMath::Sqrt(rfft[i]*rfft[i]+ifft[i]*ifft[i])/nPoints;
-    if (lmag<threshold) continue;
-    if (locMax){
-      if ( TMath::Sqrt(rfft[i-1]*rfft[i-1]+ifft[i-1]*ifft[i-1])/nPoints>lmag) continue;
-      if ( TMath::Sqrt(rfft[i+1]*rfft[i+1]+ifft[i+1]*ifft[i+1])/nPoints>lmag) continue;
-      if ( TMath::Sqrt(rfft[i-2]*rfft[i-2]+ifft[i-2]*ifft[i-2])/nPoints>lmag) continue;
-      if ( TMath::Sqrt(rfft[i+2]*rfft[i+2]+ifft[i+2]*ifft[i+2])/nPoints>lmag) continue;
-      if ( TMath::Sqrt(rfft[i-3]*rfft[i-3]+ifft[i-3]*ifft[i-3])/nPoints>lmag) continue;
-      if ( TMath::Sqrt(rfft[i+3]*rfft[i+3]+ifft[i+3]*ifft[i+3])/nPoints>lmag) continue;
-    }
-    
-    freq[current] = Float_t(i)/Float_t(nPoints);
-    //
-    re[current] = rfft[i];
-    im[current] = ifft[i];
-    mag[current]=lmag;
-    phi[current]=TMath::ATan2(ifft[i],rfft[i]);
-    current++;
-  }
-  delete [] in;
-  delete [] rfft;
-  delete [] ifft;
-  return current;
-}
-