Changes by Theodor: Speed up of clusterizer, fix of five pad cluster unfolding, faste...
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 22 Dec 2008 14:36:44 +0000 (14:36 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 22 Dec 2008 14:36:44 +0000 (14:36 +0000)
TRD/AliTRDSignalIndex.cxx
TRD/AliTRDSignalIndex.h
TRD/AliTRDarrayADC.cxx
TRD/AliTRDarraySignal.cxx
TRD/AliTRDclusterizer.cxx
TRD/AliTRDclusterizer.h

index 5251e4a..2fd9980 100644 (file)
@@ -25,7 +25,7 @@
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
-#include "TArrayI.h"
+#include "TObject.h"
 
 #include "AliLog.h"
 
@@ -40,18 +40,17 @@ AliTRDSignalIndex::AliTRDSignalIndex()
   ,fLayer(-1)
   ,fStack(-1)
   ,fSM(-1)
-  ,fIndex(NULL)
-  ,fPositionRow(0)
-  ,fPositionCol(0)
-  ,fPositionTbin(0)
-  ,fLastRow(0)
-  ,fLastCol(0)
-  ,fLastTbin(0)
+  ,fBoolIndex(NULL)
+  ,fSortedIndex(NULL)
+  ,fMaxLimit(0)
+  ,fPositionRC(0)
+  ,fSortedWasInit(kFALSE)
+  ,fCurrRow(0)
+  ,fCurrCol(0)
+  ,fCurrTbin(0)
   ,fNrows(0)
   ,fNcols(0)
   ,fNtbins(0)
-  ,fMaxLimit(0)
-  ,fResetCounters(kTRUE)
   ,fHasEntry(kFALSE)
 {
   //
@@ -69,18 +68,17 @@ AliTRDSignalIndex::AliTRDSignalIndex(Int_t nrow, Int_t ncol,Int_t ntime)
   ,fLayer(-1)
   ,fStack(-1)
   ,fSM(-1)
-  ,fIndex(NULL)
-  ,fPositionRow(0)
-  ,fPositionCol(0)
-  ,fPositionTbin(0)
-  ,fLastRow(0)
-  ,fLastCol(0)
-  ,fLastTbin(0)
+  ,fBoolIndex(NULL)
+  ,fSortedIndex(NULL)
+  ,fMaxLimit(0)
+  ,fPositionRC(0)
+  ,fSortedWasInit(kFALSE)
+  ,fCurrRow(0)
+  ,fCurrCol(0)
+  ,fCurrTbin(0)
   ,fNrows(0)
   ,fNcols(0)
   ,fNtbins(0)
-  ,fMaxLimit(0)
-  ,fResetCounters(kTRUE)
   ,fHasEntry(kFALSE)
 {
   //
@@ -98,24 +96,28 @@ AliTRDSignalIndex::AliTRDSignalIndex(const AliTRDSignalIndex &a)
   ,fLayer(a.fLayer)
   ,fStack(a.fStack)
   ,fSM(a.fSM)
-  ,fIndex(a.fIndex)
-  ,fPositionRow(a.fPositionRow)
-  ,fPositionCol(a.fPositionCol)
-  ,fPositionTbin(a.fPositionTbin)
-  ,fLastRow(a.fLastRow)
-  ,fLastCol(a.fLastCol)
-  ,fLastTbin(a.fLastTbin)
+  ,fBoolIndex(NULL)
+  ,fSortedIndex(NULL)
+  ,fMaxLimit(a.fMaxLimit)
+  ,fPositionRC(a.fPositionRC)
+  ,fSortedWasInit(a.fSortedWasInit)
+  ,fCurrRow(a.fCurrRow)
+  ,fCurrCol(a.fCurrCol)
+  ,fCurrTbin(a.fCurrTbin)
   ,fNrows(a.fNrows)
   ,fNcols(a.fNcols)
   ,fNtbins(a.fNtbins)
-  ,fMaxLimit(a.fMaxLimit)
-  ,fResetCounters(a.fResetCounters)
   ,fHasEntry(a.fHasEntry)
 {
   //
   // Copy constructor
   //
 
+  fBoolIndex = new Bool_t[fMaxLimit];
+  memcpy(fBoolIndex, a.fBoolIndex, fMaxLimit*sizeof(Bool_t));
+
+  fSortedIndex = new Short_t[2*fMaxLimit];
+  memcpy(fSortedIndex, a.fSortedIndex, 2*fMaxLimit*sizeof(Short_t));
 }
 
 //_____________________________________________________________________________
@@ -125,9 +127,14 @@ AliTRDSignalIndex::~AliTRDSignalIndex()
   // Destructor
   //
 
-  if (fIndex) {
-    delete fIndex;
-    fIndex = NULL;
+  if (fBoolIndex) {
+    delete [] fBoolIndex;
+    fBoolIndex = NULL;
+  }
+
+if (fSortedIndex) {
+    delete [] fSortedIndex;
+    fSortedIndex = NULL;
   }
 
 }
@@ -143,19 +150,31 @@ void AliTRDSignalIndex::Copy(TObject &a) const
   ((AliTRDSignalIndex &)a).fLayer         = fLayer;
   ((AliTRDSignalIndex &)a).fStack         = fStack;
   ((AliTRDSignalIndex &)a).fSM            = fSM;
-  ((AliTRDSignalIndex &)a).fIndex         = fIndex;
-  ((AliTRDSignalIndex &)a).fPositionRow   = fPositionRow;
-  ((AliTRDSignalIndex &)a).fPositionTbin  = fPositionTbin;
-  ((AliTRDSignalIndex &)a).fLastRow       = fLastRow;
-  ((AliTRDSignalIndex &)a).fLastCol       = fLastCol;
-  ((AliTRDSignalIndex &)a).fLastTbin      = fLastTbin;
+  ((AliTRDSignalIndex &)a).fMaxLimit    = fMaxLimit;
+  ((AliTRDSignalIndex &)a).fPositionRC    = fPositionRC;
+  ((AliTRDSignalIndex &)a).fSortedWasInit = fSortedWasInit;
+  ((AliTRDSignalIndex &)a).fCurrRow       = fCurrRow;
+  ((AliTRDSignalIndex &)a).fCurrCol       = fCurrCol;
+  ((AliTRDSignalIndex &)a).fCurrTbin      = fCurrTbin;
   ((AliTRDSignalIndex &)a).fNrows         = fNrows;
   ((AliTRDSignalIndex &)a).fNcols         = fNcols;
   ((AliTRDSignalIndex &)a).fNtbins        = fNtbins;
-  ((AliTRDSignalIndex &)a).fMaxLimit      = fMaxLimit;
-  ((AliTRDSignalIndex &)a).fResetCounters = fResetCounters;
   ((AliTRDSignalIndex &)a).fHasEntry      = fHasEntry;
 
+  if(((AliTRDSignalIndex &)a).fBoolIndex)
+    {
+      delete [] ((AliTRDSignalIndex &)a).fBoolIndex;
+    }
+  ((AliTRDSignalIndex &)a).fBoolIndex = new Bool_t[fMaxLimit];
+  memcpy(((AliTRDSignalIndex &)a).fBoolIndex, fBoolIndex, fMaxLimit*sizeof(Bool_t));
+
+  if(((AliTRDSignalIndex &)a).fSortedIndex)
+    {
+      delete [] ((AliTRDSignalIndex &)a).fSortedIndex;
+    }
+  ((AliTRDSignalIndex &)a).fSortedIndex = new Short_t[2*fMaxLimit];
+  memcpy(((AliTRDSignalIndex &)a).fSortedIndex, fSortedIndex, 2*fMaxLimit*sizeof(Short_t));
+
 }
 
 //_____________________________________________________________________________
@@ -171,7 +190,7 @@ AliTRDSignalIndex& AliTRDSignalIndex::operator = (const AliTRDSignalIndex& a)
 }
 
 //_____________________________________________________________________________
-void AliTRDSignalIndex::Allocate(Int_t nrow, Int_t ncol,Int_t ntime)
+void AliTRDSignalIndex::Allocate(const Int_t nrow, const Int_t ncol, const Int_t ntime)
 {
   //
   // Create the arrays
@@ -181,15 +200,22 @@ void AliTRDSignalIndex::Allocate(Int_t nrow, Int_t ncol,Int_t ntime)
   fNcols = ncol;
   fNtbins = ntime;
 
-  fMaxLimit = nrow * ncol * ntime + nrow * ncol * 2;
-  if (fIndex) {
-    delete fIndex;
-    fIndex = NULL;
+  fMaxLimit = nrow * ncol + 1;
+
+  if (fBoolIndex) {
+    delete [] fBoolIndex;
+    fBoolIndex = NULL;
+  }
+  if (fSortedIndex) {
+    delete [] fSortedIndex;
+    fSortedIndex = NULL;
   }
 
-  fIndex = new TArrayI(fMaxLimit);
-  fIndex->Reset(-1);
+  fBoolIndex = new Bool_t[fMaxLimit];
+  fSortedIndex = new Short_t[2*fMaxLimit];
 
+  ResetArrays();
   ResetCounters();
 
   fHasEntry = kFALSE;
@@ -197,6 +223,13 @@ void AliTRDSignalIndex::Allocate(Int_t nrow, Int_t ncol,Int_t ntime)
 }
 
 //_____________________________________________________________________________
+void AliTRDSignalIndex::ResetArrays()
+{
+  memset(fBoolIndex,0x00,sizeof(Bool_t)*fMaxLimit);
+  memset(fSortedIndex,0xFF,2*sizeof(Short_t)*fMaxLimit);
+}
+
+//_____________________________________________________________________________
 void AliTRDSignalIndex::Reset()
 {
   //
@@ -220,7 +253,7 @@ void AliTRDSignalIndex::ResetContent()
   // Reset the array but keep the size - no realloc
   //
 
-  fIndex->Reset(-1);
+  ResetArrays();
   ResetCounters();
 
   fHasEntry = kFALSE;
@@ -228,7 +261,7 @@ void AliTRDSignalIndex::ResetContent()
 }
 
 //_____________________________________________________________________________
-void AliTRDSignalIndex::ResetContentConditional(Int_t nrow, Int_t ncol,Int_t ntime)
+void AliTRDSignalIndex::ResetContentConditional(const Int_t nrow, const Int_t ncol, const Int_t ntime)
 {
   //
   // Reset the array but keep the size if no need to enlarge - no realloc
@@ -245,7 +278,7 @@ void AliTRDSignalIndex::ResetContentConditional(Int_t nrow, Int_t ncol,Int_t nti
     Allocate(nrow, ncol, ntime);
   }
   else {
-    fIndex->Reset(-1);
+    ResetArrays();
     ResetCounters();
     fHasEntry = kFALSE;
   }
@@ -268,11 +301,16 @@ void AliTRDSignalIndex::ClearAll()
   fNcols  = -1;
   fNtbins = -1;
 
-  if (fIndex) {
-    delete fIndex;
-    fIndex = NULL;
+  if (fBoolIndex) {
+    delete [] fBoolIndex;
+    fBoolIndex = NULL;
+  }
+
+  if (fSortedIndex) {
+    delete [] fSortedIndex;
+    fSortedIndex = NULL;
   }
-  fIndex = new TArrayI();
+  
   ResetCounters();
 
   fHasEntry = kFALSE;
@@ -280,55 +318,37 @@ void AliTRDSignalIndex::ClearAll()
 }
 
 //_____________________________________________________________________________
-void AliTRDSignalIndex::AddIndexTBin(Int_t row, Int_t col, Int_t tbin)
+void AliTRDSignalIndex::AddIndexTBin(Int_t row, Int_t col, Int_t /*tbin*/)
 {
   //
+  // This function is now obsolate, it will be deleted in future. 
+  //
   // Store the index row-column-tbin as an interesting one
   // The RC index is updated to!!!
   // This is to be used in the TRD clusterizer!
   //
+  AddIndexRC(row, col);
 
-  if (row * col * tbin + row * col * 2 >= fMaxLimit) {
-    AliError(Form("Out-of-limits fPositionCol + fNtbins %d. Limit is: %d"
-                 ,fPositionCol + fNtbins
-                 ,fMaxLimit));
-    return;
-  }
-
-  if ((row != fLastRow) || 
-      (col != fLastCol)) {
-
-    // New RC combination
-    if (fResetCounters == kTRUE) {
-      fPositionRow    = 0;
-      fPositionCol   = 1;  
-      fResetCounters = kFALSE;
-    }
-    else {
-      fPositionRow += fNtbins + 2;
-      fPositionCol += fNtbins + 2;
-    }
-
-    fPositionTbin = 1;
-
-    (*fIndex)[fPositionRow] = row;
-    (*fIndex)[fPositionCol] = col;
-    (*fIndex)[fPositionCol + fPositionTbin] = tbin;
-    ++fPositionTbin;
+}
 
-  }
-  else {
+//_____________________________________________________________________________
+void AliTRDSignalIndex::AddIndexRC(const Int_t row, const Int_t col)
+{
+  //
+  // Store the index row-column as an interesting one
+  // The RC index is updated to!!!
+  // This is to be used in the TRD clusterizer!
+  //
 
-    // Same RCT combination ?
-      
-    (*fIndex)[fPositionCol + fPositionTbin] = tbin;
-    ++fPositionTbin;      
-  
+  if (row * col + 1 >= fMaxLimit) {
+    AliError(Form("Out-of-limits row * col %d. Limit is: %d"
+                 ,row * col
+                 ,fMaxLimit));
+    return;
   }
-  
-  fLastRow  = row;
-  fLastCol  = col;
-  fLastTbin = tbin;
+  fBoolIndex[row*fNcols+col]=kTRUE;
 
   fHasEntry = kTRUE;
 
@@ -338,53 +358,44 @@ void AliTRDSignalIndex::AddIndexTBin(Int_t row, Int_t col, Int_t tbin)
 Bool_t  AliTRDSignalIndex::NextRCIndex(Int_t &row, Int_t &col)
 {
   //
-  // Return the position (index in the data array) of the next available pad
+  // Returns next used RC combination
   //
 
-  if (fResetCounters == kTRUE) {
-    fPositionRow   = 0;
-    fPositionCol   = 1;
-    fResetCounters = kFALSE;
+  if(fSortedIndex[fPositionRC]>-1){
+    row = fCurrRow = fSortedIndex[fPositionRC];
+    fPositionRC++;
+    col = fCurrCol = fSortedIndex[fPositionRC];
+    fPositionRC++;
+    return kTRUE;
   }
   else {
-    fPositionRow  += fNtbins + 2;
-    fPositionCol  += fNtbins + 2;
-  }
-
-  if (fPositionRow >= fMaxLimit) {
-    return kFALSE;
+    if(fSortedWasInit || !fHasEntry)
+      {
+        ResetCounters();
+       row = fCurrRow;
+       col = fCurrCol;
+       return kFALSE;
+      }
+    else
+      {
+       InitSortedIndex();
+       return NextRCIndex(row, col);
+      }
   }
 
-  fPositionTbin = 1;
-
-  row = (*fIndex)[fPositionRow];
-  col = (*fIndex)[fPositionCol];
-
-  if ((row > -1) && 
-      (col > -1)) { 
-    return kTRUE;
-  }
-  
-  return kFALSE;
-
 }
 
 //_____________________________________________________________________________
 Bool_t AliTRDSignalIndex::NextRCTbinIndex(Int_t &row, Int_t &col, Int_t &tbin)
 {
   //
-  // Return the position (index in the data array) of the next available tbin 
-  // within the current pad
-  //
-
-  if (fPositionRow >= fMaxLimit) {
-    return kFALSE;
-  }
+  // Returns the next tbin, or if there is no next time bin, it returns the
+  // next used RC combination.
+  //  
 
   if (NextTbinIndex(tbin)) {
-    row  = (*fIndex)[fPositionRow];
-    col  = (*fIndex)[fPositionCol];
-    fResetCounters = kFALSE;
+    row = fCurrRow;
+    col = fCurrCol;
     return kTRUE;
   }
   else {
@@ -401,24 +412,36 @@ Bool_t AliTRDSignalIndex::NextRCTbinIndex(Int_t &row, Int_t &col, Int_t &tbin)
 Bool_t AliTRDSignalIndex::NextTbinIndex(Int_t &tbin)
 {
   //
-  // Return the position (index in the data array) of the next available tbin 
-  // within the current pad
+  // Returns the next tbin of the current RC combination
   //
+  
+  if(fCurrTbin<fNtbins)
+    {
+      tbin = fCurrTbin++;
+      return kTRUE;
+    }
 
-  if ((fPositionCol + fPositionTbin >= fMaxLimit) || 
-      (fPositionTbin                >  fNtbins  )) {
-    return kFALSE;
-  }
-
-  tbin = (*fIndex)[fPositionCol + fPositionTbin];
+  return kFALSE;
 
-  if (tbin > -1) {
-    ++fPositionTbin;
-    return kTRUE;
-  }
+}
 
-  return kFALSE;
+//_____________________________________________________________________________
+void AliTRDSignalIndex::InitSortedIndex()
+{
+  //
+  // Creates the SortedIndex
+  //
 
+  fSortedWasInit = kTRUE;
+  int pos=0;
+  for(int row = 0; row < fNrows; row++)
+    for(int col = 0; col < fNcols; col++)
+      if(IsBoolIndex(row, col)){
+       fSortedIndex[pos] = row;
+       pos++;
+       fSortedIndex[pos] = col;
+       pos++;
+      }
 }
 
 //_____________________________________________________________________________
@@ -428,12 +451,8 @@ void AliTRDSignalIndex::ResetCounters()
   // Reset the counters/iterators
   //
 
-  fPositionRow   =  0;
-  fPositionCol   =  fPositionRow + 1;
-  fPositionTbin  =  1;
-  fLastRow       = -1;
-  fLastCol       = -1;
-  fLastTbin      = -1;
-  fResetCounters = kTRUE; 
-
+  fCurrRow    = -1;
+  fCurrCol    = -1;
+  fCurrTbin   = -1;
+  fPositionRC =  0;
 }
index 929c86e..5af2f40 100644 (file)
@@ -18,8 +18,6 @@
 //                                                                        //
 ////////////////////////////////////////////////////////////////////////////
 
-#include "TArrayI.h"
-
 class AliTRDSignalIndex : public TObject
 {
 
@@ -31,16 +29,19 @@ class AliTRDSignalIndex : public TObject
   virtual ~AliTRDSignalIndex();
   AliTRDSignalIndex &operator=(const AliTRDSignalIndex &d); 
 
-  virtual void     Copy(TObject &d) const; 
-  virtual void     Allocate(Int_t nrow, Int_t ncol,Int_t ntime);
+  virtual void     Copy(TObject &d) const;
+  virtual void     Allocate(const Int_t nrow, const Int_t ncol, const Int_t ntime);
 
   virtual void     Reset();
-  virtual void     ResetContentConditional(Int_t nrow, Int_t ncol,Int_t ntime);
+  virtual void     ResetContentConditional(const Int_t nrow, const Int_t ncol, const Int_t ntime);
   virtual void     ResetContent();
   virtual void     ResetCounters();
-  virtual void     ResetTbinCounter() { fPositionTbin  = 1; }
+  virtual void     ResetTbinCounter() { }
+
+          void     ResetArrays();
 
   virtual void     AddIndexTBin(Int_t row, Int_t col, Int_t tbin);
+          void     AddIndexRC(const Int_t row, const Int_t col);
 
           // Get the next pad (row and column) and return kTRUE on success
           Bool_t   NextRCIndex(Int_t &row, Int_t &col); 
@@ -49,27 +50,30 @@ class AliTRDSignalIndex : public TObject
           // Get the next active timebin and return kTRUE on success
           Bool_t   NextTbinIndex(Int_t &tbin); 
 
-          Int_t    GetCurrentRow() const  { return (*fIndex)[fPositionRow];                 }
-          Int_t    GetCurrentCol() const  { return (*fIndex)[fPositionCol];                 }
-          Int_t    GetCurrentTbin() const { return (*fIndex)[fPositionCol + fPositionTbin]; }
-  
-          // Clear the array, actually destroy and recreate w/o allocating
+          Int_t    GetCurrentRow() const  { return fCurrRow; }
+          Int_t    GetCurrentCol() const  { return fCurrCol; }
+          Int_t    GetCurrentTbin() const { return fCurrTbin; }
+
+         Bool_t   IsBoolIndex(Int_t row, Int_t col) const {return fBoolIndex[row*fNcols+col];};
+          void     InitSortedIndex();
+
+         // Clear the array, actually destroy and recreate w/o allocating
           void     ClearAll(); 
           // Return kTRUE if array allocated and there is no need to call allocate
-          Bool_t   IsAllocated() const    { if (!fIndex)                return kFALSE; 
-                                            if (fIndex->GetSize() <= 0) return kFALSE; 
-                                            else                        return kTRUE;       }
-
-          void     SetSM(Int_t ix)        { fSM      =    ix; }
-          void     SetStack(Int_t ix)     { fStack   =    ix; }
-          void     SetLayer(Int_t ix)     { fLayer   =    ix; }
-          void     SetDetNumber(Int_t ix) { fDet     =    ix; }
+          Bool_t   IsAllocated() const    { if (!fBoolIndex)    return kFALSE; 
+                                            if (fMaxLimit <= 0) return kFALSE; 
+                                            else                return kTRUE;}
+
+          void     SetSM(const Int_t ix)        { fSM      =    ix; }
+          void     SetStack(const Int_t ix)     { fStack   =    ix; }
+          void     SetLayer(const Int_t ix)     { fLayer   =    ix; }
+          void     SetDetNumber(const Int_t ix) { fDet     =    ix; }
   
   virtual Int_t    GetDetNumber() const   { return fDet;      } // Get Det number
   virtual Int_t    GetLayer() const       { return fLayer;    } // Layer position of the chamber in TRD
   virtual Int_t    GetStack() const       { return fStack;    } // Stack position of the chamber in TRD
   virtual Int_t    GetSM() const          { return fSM;       } // Super module of the TRD
-          TArrayI *GetArray() const       { return fIndex;    } // Get the tarrayi pointer for god knows what reason
+          Short_t *GetArray() const       { return fSortedIndex; } // Get the array pointer for god knows what reason
 
   virtual Bool_t   HasEntry() const       { return fHasEntry; } // Return status if has an entry
 
@@ -84,25 +88,91 @@ class AliTRDSignalIndex : public TObject
   Int_t     fStack;              //  Stack position in the full TRD
   Int_t     fSM;                 //  Super module - position in the full TRD
 
-  TArrayI  *fIndex;              //! Monitor active pads and tbins
-
-  Int_t     fPositionRow;        //  Position in the index - jumps by 1 + 1 + fNtbins
-  Int_t     fPositionCol;        //  Position in the index - jumps by 1 + 1 + fNtbins
-  Int_t     fPositionTbin;       //  Position in the tbin - goes from 0 to fNtbins
-
-  Int_t     fLastRow;            //  To keep track what is the RC combination
-  Int_t     fLastCol;            //  To keep track what is the RC combination
-  Int_t     fLastTbin;           //  To keep track what is the Tbin - will catch if raw data bogus
+  Bool_t   *fBoolIndex;          //  
+  Short_t  *fSortedIndex;        //  
+  Int_t     fMaxLimit;           //  Max number of things in the array
+  Int_t     fPositionRC;         //  Position in the SortedIndex
+  Bool_t    fSortedWasInit;      //  Was SortedIndex initialized?
 
+  Int_t     fCurrRow;            //  Last Row read out of SortedIndex
+  Int_t     fCurrCol;            //  Last Col read out of SortedIndex
+  Int_t     fCurrTbin;           //  Last outgiven Tbin
+  
   Int_t     fNrows;              //  Number of rows in the chamber
   Int_t     fNcols;              //  Number of cols in the chamber
   Int_t     fNtbins;             //  Number of tbins in the chamber
 
-  Int_t     fMaxLimit;           //  Max number of things in the array  = nrow*ncol*ntime + nrow*ncol*2
-  Bool_t    fResetCounters;      //  Reset counter status
   Bool_t    fHasEntry;           //  kTRUE flag if we have an entry 
 
   ClassDef(AliTRDSignalIndex,2)  //  Data container for one TRD detector segment
 
 };
 #endif
+
+/*
+Comment from 22 Dec 2008
+
+The structure of the Index was changed. Now no Tbin is saved anymore,
+only RC combination are saved! (reasons see below)
+
+For the readout, all tbins for a RC combination must be read out to find 
+the time bin of signal > 0.
+
+THE WRITING PROCEDURE:
+AddIndexTBin is now obsolate, use AddIndexRC instead as AddIndexTBin will
+be deleted in future.
+
+example that gives exactely the same output as before:
+as it was: 
+           AliTRDSignalIndexes *indexes;
+           AliTRDarrayADC *Signal; //or AliTRDarraySignal *Signal;
+          if(Signal->GetDataB(row, col, time)>0)
+               indexes->AddIndexTBin(row, col, time);
+
+as it should be from no on: 
+           AliTRDSignalIndexes *indexes;
+           AliTRDarrayADC *Signal; //or AliTRDarraySignal *Signal;
+          if(Signal->GetDataB(row, col, time)>0)
+               indexes->AddIndexRC(row, col);
+
+
+
+THE READING PROCEDURE:
+In most cases you can leave anything as it is.
+See more in the example.
+
+example:
+as it was: 
+           AliTRDSignalIndexes *indexes;
+           AliTRDarraySignal *Signal;
+           while(indexes->NextRCTbinIndex(row, col, time)) 
+           {...}
+
+as it should be from no on to get the exactely the same output as before: 
+           AliTRDSignalIndexes *indexes;
+           AliTRDarraySignal *Signal;
+           while(indexes->NextRCTbinIndex(row, col, time)) 
+              if(Signal->GetData(row, col, time)>0)
+                 {...}
+
+as it should be idealy:
+           AliTRDSignalIndexes *indexes;
+           AliTRDarraySignal *Signal;
+           for(time = 0; time < Ntime; time++)
+              while(indexes->NextRCIndex(row, col, time)) 
+                 if(Signal->GetData(row, col, time)>0)
+                    {...}
+
+
+REASON OF THE CHANGES:
+
+The array saved the information nicely, but it turned out that sorting 
+the array by column would have many benefits.
+I.e. it is crucial for fivePadClusters and it if much faster to allocate.
+But the sorting is not fast if the tbin is also saved.
+Moreover the tbin information was alsmost useless because, 
+whenever an RC index existed, many of the possible tbins where used.
+
+Theodor Rascanu
+
+*/
index da85ef6..91f96cc 100644 (file)
@@ -78,10 +78,7 @@ AliTRDarrayADC::AliTRDarrayADC(const AliTRDarrayADC &b)
   //
 
   fADC =  new Short_t[fNAdim];
-  for(Int_t i=0; i<fNAdim; i++)
-    {
-      fADC[i]=b.fADC[i];
-    }
+  memcpy(fADC,b.fADC, fNAdim*sizeof(Short_t));
 
 }
 
@@ -122,10 +119,8 @@ AliTRDarrayADC &AliTRDarrayADC::operator=(const AliTRDarrayADC &b)
   fNtime=b.fNtime;
   fNAdim=b.fNAdim;
   fADC = new Short_t[fNAdim];
-  for(Int_t i=0;i<fNAdim;i++)
-    {
-      fADC[i]=b.fADC[i];
-    }
+  memcpy(fADC,b.fADC, fNAdim*sizeof(Short_t));
+  
   return *this;
 
 }
@@ -149,10 +144,7 @@ void AliTRDarrayADC::Allocate(Int_t nrow, Int_t ncol, Int_t ntime)
     }
   
   fADC = new Short_t[fNAdim];
-  for(Int_t a=0; a<fNAdim; a++)
-    {
-      fADC[a]=0;
-    }
+  memset(fADC,0,sizeof(Short_t)*fNAdim);
 
 }
 
@@ -301,11 +293,8 @@ void AliTRDarrayADC::Compress()
   Int_t *longz;            
   longz = new Int_t[fNAdim];
   Int_t k=0;
-  for(Int_t i=0; i<fNAdim;i++)
-    {
-      longz[i]=0;
-      longm[i]=0;
-    }
+  memset(longz,0,sizeof(Int_t)*fNAdim);
+  memset(longm,0,sizeof(Int_t)*fNAdim);
 
   for(Int_t i=0;i<fNAdim; i++)
     {
@@ -455,11 +444,8 @@ void AliTRDarrayADC::Expand()
   longm = new Int_t[fNAdim];
   Int_t dimexp=0;
   //Initialize arrays
-  for(Int_t i=0; i<fNAdim;i++)
-    {
-      longz[i]=0;
-      longm[i]=0;
-    }
+  memset(longz,0,sizeof(Int_t)*fNAdim);
+  memset(longm,0,sizeof(Int_t)*fNAdim);
   Int_t r2=0; 
   Int_t r3=0; 
   for(Int_t i=0; i<fNAdim;i++)
index be54014..68749de 100644 (file)
@@ -26,6 +26,7 @@
 
 #include "AliTRDarraySignal.h"
 //#include "AliLog.h"
+#include "TArray.h" //for memset
 
 ClassImp(AliTRDarraySignal)
 
@@ -79,10 +80,7 @@ AliTRDarraySignal::AliTRDarraySignal(const AliTRDarraySignal &d)
   //
 
   fSignal = new Float_t[fNdim];
-  for(Int_t i=0; i<fNdim; i++)
-    {
-      fSignal[i]=d.fSignal[i];
-    }
+  memcpy(fSignal, d.fSignal, fNdim*sizeof(Float_t));
 
 }
 
@@ -123,11 +121,7 @@ AliTRDarraySignal &AliTRDarraySignal::operator=(const AliTRDarraySignal &d)
   fNtime=d.fNtime;
   fNdim=d.fNdim;
   fSignal = new Float_t[fNdim];
-
-  for(Int_t i=0; i<fNdim; i++)
-    {
-      fSignal[i]=d.fSignal[i];
-    }
+  memcpy(fSignal,d.fSignal, fNdim*sizeof(Float_t));
 
   return *this;
 
@@ -150,10 +144,7 @@ void AliTRDarraySignal::Allocate(Int_t nrow, Int_t ncol, Int_t ntime)
       delete [] fSignal;
     }
   fSignal = new Float_t[fNdim];
-  for(Int_t i=0; i<fNdim; i++)
-    {
-      fSignal[i]=0;
-    }
+  memset(fSignal,0,sizeof(Float_t)*fNdim);
 
 }
 
@@ -193,10 +184,7 @@ void AliTRDarraySignal::Compress(Float_t minval)
   Int_t k=0;
 
   //Initialize the array
-  for(Int_t i=0; i<fNdim;i++)
-    {
-      longArr[i]=0;
-    }
+  memset(longArr,0,sizeof(Int_t)*fNdim);
 
   for(Int_t i=0;i<fNdim; i++)
     {
@@ -305,10 +293,7 @@ void AliTRDarraySignal::Expand()
   Int_t *longArr; 
   longArr = new Int_t[fNdim];
   Int_t dimexp=0;
-  for(Int_t i=0; i<fNdim;i++)
-    {
-      longArr[i]=0;
-    }
+  memset(longArr,0,sizeof(Int_t)*fNdim);
 
   Int_t r2=0;
   for(Int_t i=0; i<fNdim;i++)
index 72408f8..50d0bf2 100644 (file)
@@ -71,11 +71,28 @@ AliTRDclusterizer::AliTRDclusterizer(AliTRDReconstructor *rec)
   ,fTrackletContainer(NULL)
   ,fAddLabels(kTRUE)
   ,fRawVersion(2)
-  ,fIndexesOut(NULL)
-  ,fIndexesMaxima(NULL)
   ,fTransform(new AliTRDtransform(0))
   ,fLUTbin(0)
   ,fLUT(NULL)
+  ,fDigitsIn(NULL)
+  ,fIndexes(NULL)
+  ,fADCthresh(0)
+  ,fMaxThresh(0)
+  ,fSigThresh(0)
+  ,fMinMaxCutSigma(0)
+  ,fMinLeftRightCutSigma(0)
+  ,fLayer(0)
+  ,fDet(0)
+  ,fVolid(0)
+  ,fColMax(0)
+  ,fTimeTotal(0)
+  ,fCalGainFactorROC(NULL)
+  ,fCalGainFactorDetValue(0)
+  ,fCalNoiseROC(NULL)
+  ,fCalNoiseDetValue(0)
+  ,fDigitsOut(NULL)
+  ,fClusterROC(0)
+  ,firstClusterROC(0)
 {
   //
   // AliTRDclusterizer default constructor
@@ -111,11 +128,28 @@ AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, Al
   ,fTrackletContainer(NULL)
   ,fAddLabels(kTRUE)
   ,fRawVersion(2)
-  ,fIndexesOut(NULL)
-  ,fIndexesMaxima(NULL)
   ,fTransform(new AliTRDtransform(0))
   ,fLUTbin(0)
   ,fLUT(NULL)
+  ,fDigitsIn(NULL)
+  ,fIndexes(NULL)
+  ,fADCthresh(0)
+  ,fMaxThresh(0)
+  ,fSigThresh(0)
+  ,fMinMaxCutSigma(0)
+  ,fMinLeftRightCutSigma(0)
+  ,fLayer(0)
+  ,fDet(0)
+  ,fVolid(0)
+  ,fColMax(0)
+  ,fTimeTotal(0)
+  ,fCalGainFactorROC(NULL)
+  ,fCalGainFactorDetValue(0)
+  ,fCalNoiseROC(NULL)
+  ,fCalNoiseDetValue(0)
+  ,fDigitsOut(NULL)
+  ,fClusterROC(0)
+  ,firstClusterROC(0)
 {
   //
   // AliTRDclusterizer constructor
@@ -146,11 +180,28 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
   ,fTrackletContainer(NULL)
   ,fAddLabels(kTRUE)
   ,fRawVersion(2)
-  ,fIndexesOut(NULL)
-  ,fIndexesMaxima(NULL)
   ,fTransform(NULL)
   ,fLUTbin(0)
   ,fLUT(0)
+  ,fDigitsIn(NULL)
+  ,fIndexes(NULL)
+  ,fADCthresh(0)
+  ,fMaxThresh(0)
+  ,fSigThresh(0)
+  ,fMinMaxCutSigma(0)
+  ,fMinLeftRightCutSigma(0)
+  ,fLayer(0)
+  ,fDet(0)
+  ,fVolid(0)
+  ,fColMax(0)
+  ,fTimeTotal(0)
+  ,fCalGainFactorROC(NULL)
+  ,fCalGainFactorDetValue(0)
+  ,fCalNoiseROC(NULL)
+  ,fCalNoiseDetValue(0)
+  ,fDigitsOut(NULL)
+  ,fClusterROC(0)
+  ,firstClusterROC(0)
 {
   //
   // AliTRDclusterizer copy constructor
@@ -182,16 +233,6 @@ AliTRDclusterizer::~AliTRDclusterizer()
     fTrackletContainer = NULL;
   }
 
-  if (fIndexesOut){
-    delete fIndexesOut;
-    fIndexesOut    = NULL;
-  }
-
-  if (fIndexesMaxima){
-    delete fIndexesMaxima;
-    fIndexesMaxima = NULL;
-  }
-
   if (fTransform){
     delete fTransform;
     fTransform     = NULL;
@@ -201,6 +242,11 @@ AliTRDclusterizer::~AliTRDclusterizer()
     delete [] fLUT;
     fLUT           = NULL;
   }
+  
+  if (fDigitsOut) {
+    delete fDigitsOut;
+    fDigitsOut     = NULL;
+  }
 
 }
 
@@ -234,11 +280,28 @@ void AliTRDclusterizer::Copy(TObject &c) const
   ((AliTRDclusterizer &) c).fTrackletContainer = NULL;
   ((AliTRDclusterizer &) c).fAddLabels     = fAddLabels;
   ((AliTRDclusterizer &) c).fRawVersion    = fRawVersion;
-  ((AliTRDclusterizer &) c).fIndexesOut    = NULL;
-  ((AliTRDclusterizer &) c).fIndexesMaxima = NULL;
   ((AliTRDclusterizer &) c).fTransform     = NULL;
   ((AliTRDclusterizer &) c).fLUTbin        = 0;
   ((AliTRDclusterizer &) c).fLUT           = NULL;
+  ((AliTRDclusterizer &) c).fDigitsIn      = NULL;
+  ((AliTRDclusterizer &) c).fIndexes       = NULL;
+  ((AliTRDclusterizer &) c).fADCthresh     = 0;
+  ((AliTRDclusterizer &) c).fMaxThresh     = 0;
+  ((AliTRDclusterizer &) c).fSigThresh     = 0;
+  ((AliTRDclusterizer &) c).fMinMaxCutSigma= 0;
+  ((AliTRDclusterizer &) c).fMinLeftRightCutSigma = 0;
+  ((AliTRDclusterizer &) c).fLayer         = 0;
+  ((AliTRDclusterizer &) c).fDet           = 0;
+  ((AliTRDclusterizer &) c).fVolid         = 0;
+  ((AliTRDclusterizer &) c).fColMax        = 0;
+  ((AliTRDclusterizer &) c).fTimeTotal     = 0;
+  ((AliTRDclusterizer &) c).fCalGainFactorROC = NULL;
+  ((AliTRDclusterizer &) c).fCalGainFactorDetValue = 0;
+  ((AliTRDclusterizer &) c).fCalNoiseROC   = NULL;
+  ((AliTRDclusterizer &) c).fCalNoiseDetValue = 0;
+  ((AliTRDclusterizer &) c).fDigitsOut     = NULL;
+  ((AliTRDclusterizer &) c).fClusterROC    = 0;
+  ((AliTRDclusterizer &) c).firstClusterROC= 0;
 
 }
 
@@ -444,53 +507,6 @@ Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
 }
 
 //_____________________________________________________________________________
-void AliTRDclusterizer::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
-{
-  // 
-  // Reset the helper indexes
-  //
-
-  if (fIndexesOut)
-    {
-      // carefull here - we assume that only row number may change - most probable
-      if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
-  fIndexesOut->ResetContent();
-      else
-  fIndexesOut->ResetContentConditional(indexesIn->GetNrow()
-                                          , indexesIn->GetNcol()
-                                          , indexesIn->GetNtime());
-    }
-  else
-    {
-      fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow()
-                                        , indexesIn->GetNcol() 
-                                        , indexesIn->GetNtime());
-    }
-  
-  if (fIndexesMaxima)
-    {
-      // carefull here - we assume that only row number may change - most probable
-      if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow()) 
-        {
-    fIndexesMaxima->ResetContent();
-  }
-      else
-  {
-    fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow()
-                                                , indexesIn->GetNcol()
-                                                , indexesIn->GetNtime());
-  }
-    }
-  else
-    {
-      fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow()
-                                          , indexesIn->GetNcol()
-                                          , indexesIn->GetNtime());
-    }
-
-}
-
-//_____________________________________________________________________________
 Bool_t AliTRDclusterizer::ReadDigits()
 {
   //
@@ -631,13 +647,12 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
   }
 
 
-  AliTRDrawStreamBase *pinput = AliTRDrawStreamBase::GetRawStream(rawReader);
-  AliTRDrawStreamBase &input = *pinput;
+  AliTRDrawStreamBase *input = AliTRDrawStreamBase::GetRawStream(rawReader);
 
-  AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
+  AliInfo(Form("Stream version: %s", input->IsA()->GetName()));
   
   Int_t det    = 0;
-  while ((det = input.NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
+  while ((det = input->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
     Bool_t iclusterBranch = kFALSE;
     if (fDigitsManager->GetIndexes(det)->HasEntry()){
     iclusterBranch = MakeClusters(det);
@@ -662,8 +677,8 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
   delete fDigitsManager;
   fDigitsManager = NULL;
 
-  delete pinput;
-  pinput = NULL;
+  delete input;
+  input = NULL;
 
   AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast())); 
   return kTRUE;
@@ -692,17 +707,16 @@ UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
 }
 
 //_____________________________________________________________________________
-void AliTRDclusterizer::SetPadStatus(UChar_t status, UChar_t &out){
+void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out){
   //
   // Set the pad status into out
   // First three bits are needed for the position encoding
   //
-  status = status << 3;
-  out |= status;
+  out |= status << 3;
 }
 
 //_____________________________________________________________________________
-UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding){
+UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
   //
   // return the staus encoding of the corrupted pad
   //
@@ -710,7 +724,7 @@ UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding){
 }
 
 //_____________________________________________________________________________
-Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding){
+Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
   //
   // Return the position of the corruption
   //
@@ -725,23 +739,23 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
   //
 
   // Get the digits
-  //   digits should be expanded beforehand!
+  //   digits should be expanded beforehand! 
   //   digitsIn->Expand();
-  AliTRDarrayADC *digitsIn = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod     
+  fDigitsIn = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod     
   
   // This is to take care of switched off super modules
-  if (!digitsIn->HasData()) 
+  if (!fDigitsIn->HasData()) 
     {
       return kFALSE;
     }
 
-  AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
-  if (indexesIn->IsAllocated() == kFALSE)
+  fIndexes = fDigitsManager->GetIndexes(det);
+  if (fIndexes->IsAllocated() == kFALSE)
     {
       AliError("Indexes do not exist!");
       return kFALSE;      
     }
-    
+
   AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
   if (!calibration) 
     {
@@ -749,9 +763,7 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
       return kFALSE;  
     }
 
-  // ADC thresholds
-  // There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO
-  Float_t adcThreshold   = 0; 
+  fADCthresh = 0; 
 
   if (!fReconstructor){
     AliError("Reconstructor not set\n");
@@ -760,36 +772,19 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 
   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDReconstructor::kClusterizer);
 
-  // Threshold value for the maximum
-  Float_t maxThresh            = fReconstructor->GetRecoParam()->GetClusMaxThresh();
-  // Threshold value for the digit signal
-  Float_t sigThresh            = fReconstructor->GetRecoParam()->GetClusSigThresh();
-
-  // Threshold value for the maximum ( cut noise)
-  Float_t minMaxCutSigma       = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
-  // Threshold value for the sum pad ( cut noise)
-  Float_t minLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
-
-  // Iteration limit for unfolding procedure
-  const Float_t kEpsilon = 0.01;             
-  const Int_t   kNclus   = 3;  
-  const Int_t   kNsig    = 5;
+  fMaxThresh            = fReconstructor->GetRecoParam()->GetClusMaxThresh();
+  fSigThresh            = fReconstructor->GetRecoParam()->GetClusSigThresh();
+  fMinMaxCutSigma       = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
+  fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
 
-  Int_t    iUnfold       = 0;  
-  Double_t ratioLeft     = 1.0;
-  Double_t ratioRight    = 1.0;
-
-  Double_t padSignal[kNsig];   
-  Double_t clusterSignal[kNclus];
-
-  Int_t istack  = indexesIn->GetStack();
-  Int_t ilayer  = indexesIn->GetLayer();
-  Int_t isector = indexesIn->GetSM();
+  Int_t istack  = fIndexes->GetStack();
+  fLayer  = fIndexes->GetLayer();
+  Int_t isector = fIndexes->GetSM();
 
   // Start clustering in the chamber
 
-  Int_t idet  = AliTRDgeometry::GetDetector(ilayer,istack,isector);
-  if (idet != det) {
+  fDet  = AliTRDgeometry::GetDetector(fLayer,istack,isector);
+  if (fDet != det) {
     AliError("Strange Detector number Missmatch!");
     return kFALSE;
   }
@@ -797,329 +792,327 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
   // TRD space point transformation
   fTransform->SetDetector(det);
 
-  Int_t    iGeoLayer  = AliGeomManager::kTRD1 + ilayer;
+  Int_t    iGeoLayer  = AliGeomManager::kTRD1 + fLayer;
   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
-  UShort_t volid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule); 
+  fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule); 
 
-  Int_t nColMax    = digitsIn->GetNcol();
-  Int_t nRowMax    = digitsIn->GetNrow();
-  Int_t nTimeTotal = digitsIn->GetNtime();
+  fColMax    = fDigitsIn->GetNcol();
+  Int_t nRowMax    = fDigitsIn->GetNrow();
+  fTimeTotal = fDigitsIn->GetNtime();
 
   // Detector wise calibration object for the gain factors
-  const AliTRDCalDet           *calGainFactorDet      = calibration->GetGainFactorDet();
+  const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
   // Calibration object with pad wise values for the gain factors
-  AliTRDCalROC                 *calGainFactorROC      = calibration->GetGainFactorROC(idet);
+  fCalGainFactorROC      = calibration->GetGainFactorROC(fDet);
   // Calibration value for chamber wise gain factor
-  Float_t                       calGainFactorDetValue = calGainFactorDet->GetValue(idet);
+  fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
 
   // Detector wise calibration object for the noise
-  const AliTRDCalDet           *calNoiseDet           = calibration->GetNoiseDet();
+  const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
   // Calibration object with pad wise values for the noise
-  AliTRDCalROC                 *calNoiseROC           = calibration->GetNoiseROC(idet);
+  fCalNoiseROC           = calibration->GetNoiseROC(fDet);
   // Calibration value for chamber wise noise
-  Float_t                       calNoiseDetValue      = calNoiseDet->GetValue(idet);
-
-  Int_t nClusters = 0;
-
-  AliTRDarraySignal *digitsOut = new AliTRDarraySignal(nRowMax, nColMax, nTimeTotal);
-  AliTRDarrayADC     padStatus(nRowMax, nColMax, nTimeTotal);
+  fCalNoiseDetValue      = calNoiseDet->GetValue(fDet);
 
-  ResetHelperIndexes(indexesIn);
+  if(fDigitsOut) delete fDigitsOut;
+  fDigitsOut = new AliTRDarraySignal(nRowMax, fColMax, fTimeTotal);
+  
+  firstClusterROC = -1;
+  fClusterROC     =  0;
 
   // Apply the gain and the tail cancelation via digital filter
-  TailCancelation(digitsIn
-                 ,digitsOut  
-                 ,indexesIn
-                 ,fIndexesOut
-                 ,nTimeTotal
-                 ,adcThreshold
-                 ,calGainFactorROC
-                 ,calGainFactorDetValue);      
-  
-  Int_t row  = 0;
-  Int_t col  = 0;
-  Int_t time = 0;
-  Int_t iPad = 0;
-    
-  UChar_t status[3]={0, 0, 0}, ipos = 0;
-  fIndexesOut->ResetCounters();
-  Int_t nMaximas = 0, nCorrupted = 0;
-  while (fIndexesOut->NextRCTbinIndex(row, col, time)) {
-    // reset pad status
-    ipos = 0; for(Int_t is=3; is--;) status[is] = 0;
+  TailCancelation();   
 
-    Float_t signalM = TMath::Abs(digitsOut->GetData(row,col,time));
-    status[1] = digitsIn->GetPadStatus(row,col,time);
-    if(status[1]) SETBIT(ipos, AliTRDcluster::kMaskedCenter);
+  ClusterizerStruct curr, last;
+  last.Row = -1;
+  Int_t nMaximas = 0, nCorrupted = 0;
+  Double_t Ratio = 1;
 
-    if(signalM < maxThresh) continue; 
+  // Here the clusterfining is happening
+  
+  for(curr.Time = 0; curr.Time < fTimeTotal; curr.Time++)
+    while(fIndexes->NextRCIndex(curr.Row, curr.Col))
+      if(IsMaximum(curr.Row, curr.Col, curr.Time, curr.padStatus, &curr.Signals[0]))
+       {
+         if(last.Row>-1)
+               {
+                 last.Signals[0] *= Ratio;
+                 if(curr.Row==last.Row && curr.Col==last.Col+2)
+                   {
+                     if(IsFivePadCluster(last.Row, last.Col, last.Time, &last.Signals[0], &curr.Signals[0], Ratio))
+                       {
+                         last.Signals[2] *= Ratio;
+                         Ratio = 1 - Ratio;
+                       }else Ratio = 1;
+                   }else Ratio = 1;
+                 CreateCluster(last.Row, last.Col, last.Time, &last.Signals[0], last.padStatus);
+               }
+         last=curr;
+       }
+  if(last.Row>-1)
+    {
+      last.Signals[0] *= Ratio;
+      CreateCluster(last.Row, last.Col, last.Time, &last.Signals[0], last.padStatus);
+    }
 
-    Float_t  noiseMiddleThresh = minMaxCutSigma*calNoiseDetValue*calNoiseROC->GetValue(col,row);
-    if (signalM < noiseMiddleThresh) continue;
+  if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 2){
+    (*fDebugStream) << "MakeClusters"
+                   << "Detector="   << det
+                   << "NMaxima="    << nMaximas
+                   << "NClusters="  << fClusterROC
+                   << "NCorrupted=" << nCorrupted
+                   << "\n";
+  }
 
-    if (col + 1 >= nColMax || col-1 < 0) continue;
-    
-    Float_t signalL = TMath::Abs(digitsOut->GetData(row,col+1,time));
-    status[0] = digitsIn->GetPadStatus(row,col+1,time);
-    if(status[0]) SETBIT(ipos, AliTRDcluster::kMaskedLeft);
-    
-    Float_t signalR = TMath::Abs(digitsOut->GetData(row,col-1,time));
-    status[2] = digitsIn->GetPadStatus(row,col-1,time);
-    if(status[2]) SETBIT(ipos, AliTRDcluster::kMaskedRight);
-    
-    // reject candidates with more than 1 problematic pad
-    if(ipos == 3 || ipos > 4) continue;
-    
-    if (!status[1]) { // good central pad
-      if (!ipos) { // all pads are OK
-        if ((signalL <= signalM) && (signalR <  signalM)) {
-          if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
-             Float_t  noiseSumThresh = minLeftRightCutSigma
-                                     * calNoiseDetValue
-                                     * calNoiseROC->GetValue(col,row);
-            if ((signalL+signalR+signalM) >= noiseSumThresh) {
-              // Maximum found, mark the position by a negative signal
-              digitsOut->SetData(row,col,time,-signalM);
-              fIndexesMaxima->AddIndexTBin(row,col,time);
-              padStatus.SetData(row, col, time, ipos); // No corruption
-            }
-          }
-        }
-      } else { // one of the neighbouring pads are bad
-        if (status[0] && signalR < signalM && signalR >= sigThresh) {
-          digitsOut->SetData(row,col,time,-signalM);
-          digitsOut->SetData(row, col, time+1, 0.);
-          fIndexesMaxima->AddIndexTBin(row,col,time);
-          SetPadStatus(status[0], ipos);
-          padStatus.SetData(row, col, time, ipos);
-        } 
-        else if (status[2] && signalL <= signalM && signalL >= sigThresh) {
-          digitsOut->SetData(row,col,time,-signalM);
-          digitsOut->SetData(row, col, time-1, 0.);
-          fIndexesMaxima->AddIndexTBin(row,col,time);
-          SetPadStatus(status[2], ipos);
-          padStatus.SetData(row, col, time, ipos);
-        }
-      }
-    } 
-    else { // wrong maximum pad
-      if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
-        // Maximum found, mark the position by a negative signal
-        digitsOut->SetData(row,col,time,-maxThresh);
-        fIndexesMaxima->AddIndexTBin(row,col,time);
-        SetPadStatus(status[1], ipos);
-        padStatus.SetData(row, col, time, ipos);
-      }
-    }
+  if (fAddLabels) {
+    AddLabels(fDet,firstClusterROC,fClusterROC);
   }
 
-  // The index to the first cluster of a given ROC
-  Int_t firstClusterROC = -1;
-  // The number of cluster in a given ROC
-  Int_t nClusterROC     =  0;
+  return kTRUE;
 
-  // Now check the maxima and calculate the cluster position
-  fIndexesMaxima->ResetCounters();
-  while (fIndexesMaxima->NextRCTbinIndex(row, col, time)) {
+}
 
-    // Maximum found ?             
-    if (digitsOut->GetData(row,col,time) < 0.0) {
+//_____________________________________________________________________________
+Bool_t AliTRDclusterizer::IsMaximum(const Int_t row, const Int_t col, const Int_t time, 
+                                   UChar_t &pasStatus, Double_t *const Signals) 
+{
+  //
+  // Returns true if this row,col,time combination is a maximum. 
+  // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
+  //
 
-      for (iPad = 0; iPad < kNclus; iPad++) {
-        Int_t iPadCol = col - 1 + iPad;
-        clusterSignal[iPad] = TMath::Abs(digitsOut->GetData(row,iPadCol,time));
-      }
+  Signals[1] = fDigitsOut->GetData(row,col,time);
+  if(Signals[1] < fMaxThresh) return kFALSE;
 
-      // Count the number of pads in the cluster
-      Int_t nPadCount = 0;
-      Int_t ii;
-      // Look to the right
-      ii = 0;
-      while (TMath::Abs(digitsOut->GetData(row,col-ii  ,time)) >= sigThresh) {
-        nPadCount++;
-        ii++;
-        if (col-ii   <        0) break;
-      }
-      // Look to the left
-      ii = 0;
-      while (TMath::Abs(digitsOut->GetData(row,col+ii+1,time)) >= sigThresh) {
-        nPadCount++;
-        ii++;
-        if (col+ii+1 >= nColMax) break;
-      }
-      nClusters++;
+  Float_t  noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(col,row);
+  if (Signals[1] < noiseMiddleThresh) return kFALSE;
 
-      // Look for 5 pad cluster with minimum in the middle
-      Bool_t fivePadCluster = kFALSE;
-      if (col < (nColMax - 3)){
-        if (digitsOut->GetData(row,col+2,time) < 0) {
-          fivePadCluster = kTRUE;
-        }
-        if ((fivePadCluster) && (col < (nColMax - 5))) {
-          if (digitsOut->GetData(row,col+4,time) >= sigThresh) {
-            fivePadCluster = kFALSE;
-          }
-        }
-        if ((fivePadCluster) && (col >             1)) {
-          if (digitsOut->GetData(row,col-2,time) >= sigThresh) {
-            fivePadCluster = kFALSE;
-          }
-        }
-      }
+  if (col + 1 >= fColMax || col < 1) return kFALSE;
+  UChar_t status[3]={0, 0, 0};
+  pasStatus = 0;
 
-      // 5 pad cluster
-      // Modify the signal of the overlapping pad for the left part 
-      // of the cluster which remains from a previous unfolding
-      if (iUnfold) {
-        clusterSignal[0] *= ratioLeft;
-        iUnfold = 0;
-      }
+  status[1] = fDigitsIn->GetPadStatus(row,col,time);
+  //if(status[1]) SETBIT(pasStatus, AliTRDcluster::kMaskedCenter);//TR: mod: this is already done by SetPadStatus
 
-      // Unfold the 5 pad cluster
-      if (fivePadCluster) {
-        for (iPad = 0; iPad < kNsig; iPad++) {
-          padSignal[iPad] = TMath::Abs(digitsOut->GetData(row
-                                                         ,col-1+iPad
-                                                         ,time));
-        }
-        // Unfold the two maxima and set the signal on 
-        // the overlapping pad to the ratio
-        ratioRight        = Unfold(kEpsilon,ilayer,padSignal);
-        ratioLeft         = 1.0 - ratioRight; 
-        clusterSignal[2] *= ratioRight;
-        iUnfold = 1;
+  Signals[2] = fDigitsOut->GetData(row,col+1,time);
+  status[2] = fDigitsIn->GetPadStatus(row,col+1,time);
+  //if(status[2]) SETBIT(pasStatus, AliTRDcluster::kMaskedLeft);//TR: mod: this is already done by SetPadStatus
+    
+  Signals[0] = fDigitsOut->GetData(row,col-1,time);
+  status[0] = fDigitsIn->GetPadStatus(row,col-1,time);
+  //if(status[0]) SETBIT(pasStatus, AliTRDcluster::kMaskedRight);//TR: mod: this is already done by SetPadStatus
+    
+  // reject candidates with more than 1 problematic pad
+  if(pasStatus >= 3) return kFALSE;
+    
+  if (!status[1]) { // good central pad
+    if (!pasStatus) { // all pads are OK
+      if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
+       if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
+         Float_t  noiseSumThresh = fMinLeftRightCutSigma
+           * fCalNoiseDetValue
+           * fCalNoiseROC->GetValue(col,row);
+         if ((Signals[2]+Signals[0]+Signals[1]) >= noiseSumThresh)
+           return kTRUE;
+       }
       }
-
-      // The position of the cluster in COL direction relative to the center pad (pad units)
-      Double_t clusterPosCol = 0.0;
-      if (fReconstructor->GetRecoParam()->IsLUT()) {
-        // Calculate the position of the cluster by using the
-        // lookup table method
-        clusterPosCol = LUTposition(ilayer,clusterSignal[0]
-                                          ,clusterSignal[1]
-                                          ,clusterSignal[2]);
+    } else { // one of the neighbouring pads are bad
+      if (status[2] && Signals[0] < Signals[1] && Signals[0] >= fSigThresh) { 
+       fDigitsOut->SetData(row, col+1, time, 0.);//TR: mod: was: SetData(row, col, time+1, 0.)
+       SetPadStatus(status[2], pasStatus);
+       return kTRUE;
       } 
-      else {
-        // Calculate the position of the cluster by using the
-        // center of gravity method
-        for (Int_t i = 0; i < kNsig; i++) {
-          padSignal[i] = 0.0;
-        }
-        padSignal[2] = TMath::Abs(digitsOut->GetData(row,col  ,time)); // Central pad
-        padSignal[3] = TMath::Abs(digitsOut->GetData(row,col+1,time)); // Left    pad
-        padSignal[1] = TMath::Abs(digitsOut->GetData(row,col-1,time)); // Right   pad
-        if ((col >           2) && 
-            (TMath::Abs(digitsOut->GetData(row,col-2,time)) < padSignal[1])) {
-         padSignal[4] = TMath::Abs(digitsOut->GetData(row,col-2,time));
-        }
-        if ((col < nColMax - 3) &&
-            (TMath::Abs(digitsOut->GetData(row,col+2,time)) < padSignal[3])) {
-          padSignal[0] = TMath::Abs(digitsOut->GetData(row,col+2,time));
-        }
-        clusterPosCol = GetCOG(padSignal);
+      else if (status[0] && Signals[2] <= Signals[1] && Signals[2] >= fSigThresh) { 
+       fDigitsOut->SetData(row, col-1, time, 0.);//TR: mod: was: SetData(row, col, time-1, 0.)
+       SetPadStatus(status[0], pasStatus);
+       return kTRUE;
       }
+    }
+  } 
+  else { // wrong maximum pad
+    if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
+      fDigitsOut->SetData(row,col,time,fMaxThresh);
+      SetPadStatus(status[1], pasStatus);
+      return kTRUE;
+    }
+  }
+  return kFALSE;
+}
 
-      // Store the amplitudes of the pads in the cluster for later analysis
-      // and check whether one of these pads is masked in the database
-      Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
-      for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
-        if ((jPad <          0) || 
-            (jPad >= nColMax-1)) {
-          continue;
-        }
-        signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetData(row,jPad,time)));
-      }
-
-      // Transform the local cluster coordinates into calibrated 
-      // space point positions defined in the local tracking system.
-      // Here the calibration for T0, Vdrift and ExB is applied as well.
-      Double_t clusterXYZ[6];
-      clusterXYZ[0] = clusterPosCol;
-      clusterXYZ[1] = clusterSignal[2];
-      clusterXYZ[2] = clusterSignal[1];
-      clusterXYZ[3] = clusterSignal[0];
-      clusterXYZ[4] = 0.0;
-      clusterXYZ[5] = 0.0;
-      Int_t    clusterRCT[3];
-      clusterRCT[0] = row;
-      clusterRCT[1] = col;
-      clusterRCT[2] = 0;
-
-      Bool_t out = kTRUE;
-      if (fTransform->Transform(clusterXYZ,clusterRCT,((UInt_t) time),out,0)) {
-
-        // Add the cluster to the output array
-        // The track indices will be stored later 
-        Float_t clusterPos[3];
-        clusterPos[0] = clusterXYZ[0];
-        clusterPos[1] = clusterXYZ[1];
-        clusterPos[2] = clusterXYZ[2];
-        Float_t clusterSig[2];
-        clusterSig[0] = clusterXYZ[4];
-        clusterSig[1] = clusterXYZ[5];
-        Double_t clusterCharge  = clusterXYZ[3];
-        Char_t   clusterTimeBin = ((Char_t) clusterRCT[2]);
-
-        Int_t n = RecPoints()->GetEntriesFast();
-        AliTRDcluster *cluster = new((*RecPoints())[n]) AliTRDcluster(
-          idet,
-          clusterCharge, clusterPos, clusterSig,
-          0x0,
-          ((Char_t) nPadCount),
-          signals,
-          ((UChar_t) col), ((UChar_t) row), ((UChar_t) time),
-          clusterTimeBin, clusterPosCol,
-          volid);
-        cluster->SetInChamber(!out);
-
-        UChar_t maskPosition = GetCorruption(padStatus.GetData(row, col, time));
-        UChar_t padstatus = GetPadStatus(padStatus.GetData(row, col, time));
-        if (maskPosition) { 
-          cluster->SetPadMaskedPosition(maskPosition);
-          cluster->SetPadMaskedStatus(padstatus);
-        }
-
-        // Temporarily store the row, column and time bin of the center pad
-        // Used to later on assign the track indices
-        cluster->SetLabel( row,0);
-        cluster->SetLabel( col,1);
-        cluster->SetLabel(time,2);
+//_____________________________________________________________________________
+Bool_t AliTRDclusterizer::IsFivePadCluster(const Int_t row, const Int_t col, const Int_t time, 
+                                          Double_t *SignalsThisMax, Double_t *SignalsNeighbourMax, Double_t &ratio)
+{
+  //
+  // Look for 5 pad cluster with minimum in the middle
+  // Gives back the ratio
+  //
 
-        // Store the index of the first cluster in the current ROC
-        if (firstClusterROC < 0) {
-          firstClusterROC = RecPoints()->GetEntriesFast() - 1;
-        }
+  if (col < fColMax - 3){
+    if (col < fColMax - 5){
+      if (fDigitsOut->GetData(row,col+4,time) >= fSigThresh)
+       return kFALSE;
+    }
+    if (col > 1) {
+      if (fDigitsOut->GetData(row,col-2,time) >= fSigThresh)
+       return kFALSE;
+       }
+    
+  //if (fSignalsThisMax[1] >= 0){ //TR: mod
 
-        // Count the number of cluster in the current ROC
-        nClusterROC++;
+    const Float_t kEpsilon = 0.01;
+    Double_t padSignal[5] = {SignalsThisMax[0],SignalsThisMax[1],SignalsThisMax[2],
+                            SignalsNeighbourMax[1], SignalsNeighbourMax[2]};
+    
+    // Unfold the two maxima and set the signal on 
+    // the overlapping pad to the ratio
+    ratio = Unfold(kEpsilon,fLayer,padSignal);
+    return kTRUE;
+  }
+  return kFALSE;
+}
 
-      } // if: Transform ok ?
+//_____________________________________________________________________________
+void AliTRDclusterizer::CreateCluster(const Int_t row, const Int_t col, const Int_t time, 
+                                       const Double_t* const clusterSignal, const UChar_t pasStatus)
+{
+  //
+  // Creates a cluster at the given position and saves it in fRecPoint
+  //
 
-    } // if: Maximum found ?
+  const Int_t kNsig = 5;
+  Double_t padSignal[kNsig];
 
+  // The position of the cluster in COL direction relative to the center pad (pad units)
+  Double_t clusterPosCol = 0.0;
+  if (fReconstructor->GetRecoParam()->IsLUT()) {
+    // Calculate the position of the cluster by using the
+    // lookup table method
+    clusterPosCol = LUTposition(fLayer,clusterSignal[0]
+                               ,clusterSignal[1]
+                               ,clusterSignal[2]);
+  } 
+  else {
+    // Calculate the position of the cluster by using the
+    // center of gravity method
+    padSignal[1] = clusterSignal[0];
+    padSignal[2] = clusterSignal[1];
+    padSignal[3] = clusterSignal[2];
+    if(col > 2){
+      padSignal[0] = fDigitsOut->GetData(row,col-2,time);
+      if(padSignal[0]>= padSignal[1])
+       padSignal[0] = 0;
+    }
+    if(col < fColMax - 3){
+      padSignal[4] = fDigitsOut->GetData(row,col+2,time);
+      if(padSignal[4]>= padSignal[3])
+       padSignal[4] = 0;
+    }
+    clusterPosCol = GetCOG(padSignal);
   }
 
-  if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 2){
-    (*fDebugStream) << "MakeClusters"
-                   << "Detector="   << det
-                   << "NMaxima="    << nMaximas
-                   << "NClusters="  << nClusterROC
-                   << "NCorrupted=" << nCorrupted
-                   << "\n";
+  // Count the number of pads in the cluster
+  Int_t nPadCount = 1;
+  // Look to the right
+  Int_t ii = 1;
+  while (fDigitsOut->GetData(row, col-ii, time) >= fSigThresh) {
+    nPadCount++;
+    ii++;
+    if (col-ii < 0) break;
+  }
+  // Look to the left
+  ii = 1;
+  while (fDigitsOut->GetData(row, col+ii, time) >= fSigThresh) {
+    nPadCount++;
+    ii++;
+    if (col+ii >= fColMax) break;
   }
 
-  delete digitsOut;
+  // Store the amplitudes of the pads in the cluster for later analysis
+  // and check whether one of these pads is masked in the database
+  Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
+  for(Int_t i = 0; i++; i<3)
+    signals[i+2] = TMath::Nint(clusterSignal[i]);
+  for(Int_t i = 0; i++; i<2)
+    {
+      if(col+i >= 3)
+       signals[i] = TMath::Nint(fDigitsOut->GetData(row,col-3+i,time));
+      if(col+3-i < fColMax)
+       signals[7-i] = TMath::Nint(fDigitsOut->GetData(row,col+3-i,time));
+    }
+  /*for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
+    if ((jPad >= 0) && (jPad < fColMax))
+      signals[jPad-col+3] = TMath::Nint(fDigitsOut->GetData(row,jPad,time));
+      }*/
+
+  // Transform the local cluster coordinates into calibrated 
+  // space point positions defined in the local tracking system.
+  // Here the calibration for T0, Vdrift and ExB is applied as well.
+  Double_t clusterXYZ[6];
+  clusterXYZ[0] = clusterPosCol;
+  clusterXYZ[1] = clusterSignal[2];
+  clusterXYZ[2] = clusterSignal[1];
+  clusterXYZ[3] = clusterSignal[0];
+  clusterXYZ[4] = 0.0;
+  clusterXYZ[5] = 0.0;
+  Int_t    clusterRCT[3];
+  clusterRCT[0] = row;
+  clusterRCT[1] = col;
+  clusterRCT[2] = 0;
+
+  Bool_t out = kTRUE;
+  if (fTransform->Transform(clusterXYZ,clusterRCT,((UInt_t) time),out,0)) {
+
+    // Add the cluster to the output array
+    // The track indices will be stored later 
+    Float_t clusterPos[3];
+    clusterPos[0] = clusterXYZ[0];
+    clusterPos[1] = clusterXYZ[1];
+    clusterPos[2] = clusterXYZ[2];
+    Float_t clusterSig[2];
+    clusterSig[0] = clusterXYZ[4];
+    clusterSig[1] = clusterXYZ[5];
+    Double_t clusterCharge  = clusterXYZ[3];
+    Char_t   clusterTimeBin = ((Char_t) clusterRCT[2]);
+
+    Int_t n = RecPoints()->GetEntriesFast();
+    AliTRDcluster *cluster = new((*RecPoints())[n]) AliTRDcluster(
+                                                                 fDet,
+                                                                 clusterCharge, clusterPos, clusterSig,
+                                                                 0x0,
+                                                                 ((Char_t) nPadCount),
+                                                                 signals,
+                                                                 ((UChar_t) col), ((UChar_t) row), ((UChar_t) time),
+                                                                 clusterTimeBin, clusterPosCol,
+                                                                 fVolid);
+    cluster->SetInChamber(!out);
+
+    UChar_t maskPosition = GetCorruption(pasStatus);
+    UChar_t padstatus = GetPadStatus(pasStatus);
+    if (maskPosition) { 
+      cluster->SetPadMaskedPosition(maskPosition);
+      cluster->SetPadMaskedStatus(padstatus);
+    }
+
+    // Temporarily store the row, column and time bin of the center pad
+    // Used to later on assign the track indices
+    cluster->SetLabel(row, 0);
+    cluster->SetLabel(col, 1);
+    cluster->SetLabel(time,2);
 
-  if (fAddLabels) {
-    AddLabels(idet,firstClusterROC,nClusterROC);
-  }
+    // Store the index of the first cluster in the current ROC
+    if (firstClusterROC < 0) {
+      firstClusterROC = RecPoints()->GetEntriesFast() - 1;
+    }
 
-  return kTRUE;
+    // Count the number of cluster in the current ROC
+    fClusterROC++;
+  }
 
 }
 
 //_____________________________________________________________________________
-Bool_t AliTRDclusterizer::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
+Bool_t AliTRDclusterizer::AddLabels(const Int_t idet, const Int_t firstClusterROC, const Int_t nClusterROC)
 {
   //
   // Add the track indices to the found clusters
@@ -1211,7 +1204,7 @@ Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
 }
 
 //_____________________________________________________________________________
-Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSignal)
+Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSignal) const
 {
   //
   // Method to unfold neighbouring maxima.
@@ -1244,14 +1237,14 @@ Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSigna
 
     // Cluster position according to charge ratio
     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) 
-                      / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
+                      / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) 
                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
 
     // Set cluster charge ratio
-    irc = calibration->PadResponse(1.0,maxLeft ,layer,newSignal);
+    irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
     Double_t ampLeft  = padSignal[1] / newSignal[1];
-    irc = calibration->PadResponse(1.0,maxRight,layer,newSignal);
+    irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
     Double_t ampRight = padSignal[3] / newSignal[1];
 
     // Apply pad response to parameters
@@ -1269,40 +1262,33 @@ Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSigna
 }
 
 //_____________________________________________________________________________
-void AliTRDclusterizer::TailCancelation(AliTRDarrayADC *digitsIn
-                                     , AliTRDarraySignal *digitsOut   
-                                     , AliTRDSignalIndex *indexesIn
-                                     , AliTRDSignalIndex *indexesOut
-                                     , Int_t nTimeTotal
-                                     , Float_t adcThreshold
-                                     , AliTRDCalROC *calGainFactorROC
-                                     , Float_t calGainFactorDetValue)
+void AliTRDclusterizer::TailCancelation()
 {
   //
   // Applies the tail cancelation and gain factors: 
-  // Transform digitsIn to digitsOut
+  // Transform fDigitsIn to fDigitsOut
   //
 
   Int_t iRow  = 0;
   Int_t iCol  = 0;
   Int_t iTime = 0;
 
-  Double_t *inADC  = new Double_t[nTimeTotal];  // ADC data before tail cancellation
-  Double_t *outADC = new Double_t[nTimeTotal];  // ADC data after tail cancellation
-  indexesIn->ResetCounters();
+  Double_t *inADC  = new Double_t[fTimeTotal];  // ADC data before tail cancellation
+  Double_t *outADC = new Double_t[fTimeTotal];  // ADC data after tail cancellation
+  fIndexes->ResetCounters();
   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDReconstructor::kClusterizer);
-  while (indexesIn->NextRCIndex(iRow, iCol))
+  while(fIndexes->NextRCIndex(iRow, iCol))
     {
-      Float_t  calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
-      Double_t gain                  = calGainFactorDetValue 
-                                    * calGainFactorROCValue;
+      Float_t  fCalGainFactorROCValue = fCalGainFactorROC->GetValue(iCol,iRow);
+      Double_t gain                  = fCalGainFactorDetValue 
+                                    * fCalGainFactorROCValue;
 
       Bool_t corrupted = kFALSE;
-      for (iTime = 0; iTime < nTimeTotal; iTime++) 
+      for (iTime = 0; iTime < fTimeTotal; iTime++) 
         {        
           // Apply gain gain factor
-          inADC[iTime]   = digitsIn->GetDataB(iRow,iCol,iTime);
-          if (digitsIn->GetPadStatus(iRow, iCol, iTime)) corrupted = kTRUE;
+          inADC[iTime]   = fDigitsIn->GetDataB(iRow,iCol,iTime);
+          if (fDigitsIn->GetPadStatus(iRow, iCol, iTime)) corrupted = kTRUE;
           inADC[iTime]  /= gain;
           outADC[iTime]  = inADC[iTime];
          if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 7){
@@ -1321,22 +1307,15 @@ void AliTRDclusterizer::TailCancelation(AliTRDarrayADC *digitsIn
         {
           // Apply the tail cancelation via the digital filter
           // (only for non-coorupted pads)
-          if (fReconstructor->GetRecoParam() ->IsTailCancelation()) 
-            {
-              DeConvExp(inADC,outADC,nTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
-            } 
+          if (fReconstructor->GetRecoParam()->IsTailCancelation()) 
+           DeConvExp(inADC,outADC,fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
         }
 
-      indexesIn->ResetTbinCounter();
-
-      while (indexesIn->NextTbinIndex(iTime))
+      for(iTime = 0; iTime < fTimeTotal; iTime++)//while (fIndexes->NextTbinIndex(iTime))
         {
           // Store the amplitude of the digit if above threshold
-          if (outADC[iTime] > adcThreshold) 
-            {
-              digitsOut->SetData(iRow,iCol,iTime,outADC[iTime]);
-              indexesOut->AddIndexTBin(iRow,iCol,iTime);
-            }    
+          if (outADC[iTime] > fADCthresh) 
+           fDigitsOut->SetData(iRow,iCol,iTime,outADC[iTime]);
         } // while itime
 
     } // while irow icol
@@ -1349,8 +1328,8 @@ void AliTRDclusterizer::TailCancelation(AliTRDarrayADC *digitsIn
 }
 
 //_____________________________________________________________________________
-void AliTRDclusterizer::DeConvExp(Double_t *source, Double_t *target
-                                , Int_t n, Int_t nexp) 
+void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
+                                 ,const Int_t n, const Int_t nexp) 
 {
   //
   // Tail cancellation by deconvolution for PASA v4 TRF
index 80a5869..b4b7f2b 100644 (file)
@@ -69,10 +69,9 @@ class AliTRDclusterizer : public TNamed
   virtual Bool_t   MakeClusters();
   virtual Bool_t   MakeClusters(Int_t det);
 
-  virtual Bool_t   AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC);
-  virtual Bool_t   SetAddLabels(Bool_t kset) { fAddLabels = kset; 
-            return fAddLabels;  } // should we assign labels to clusters
-  virtual void     SetRawVersion(Int_t iver) { fRawVersion = iver; } // set the expected raw data version
+  virtual Bool_t   AddLabels(const Int_t idet, const Int_t firstClusterROC, const Int_t nClusterROC);
+  virtual Bool_t   SetAddLabels(const Bool_t kset) { fAddLabels = kset; return fAddLabels;  } // should we assign labels to clusters
+  virtual void     SetRawVersion(const Int_t iver) { fRawVersion = iver; } // set the expected raw data version
   void             SetReconstructor(const AliTRDReconstructor *rec) {fReconstructor = rec;}
   static UChar_t   GetStatus(Short_t &signal);
 
@@ -81,25 +80,25 @@ class AliTRDclusterizer : public TNamed
 
  protected:
 
-  void             DeConvExp(Double_t *source, Double_t *target
-                           , Int_t nTimeTotal, Int_t nexp);
-  void             TailCancelation(AliTRDarrayADC *digitsIn
-                                , AliTRDarraySignal *digitsOut 
-                                , AliTRDSignalIndex *indexesIn
-                                , AliTRDSignalIndex *indexesOut
-                                , Int_t nTimeTotal
-                                , Float_t ADCthreshold
-                                , AliTRDCalROC *calGainFactorROC
-                                , Float_t calGainFactorDetValue);
-  virtual Double_t Unfold(Double_t eps, Int_t layer, Double_t *padSignal);
+  void             DeConvExp (const Double_t *const source, Double_t *const target
+                            ,const Int_t nTimeTotal, const Int_t nexp);
+  void             TailCancelation();
+
+  virtual Double_t Unfold(Double_t eps, Int_t layer, Double_t *padSignal) const;
           Double_t GetCOG(Double_t signal[5]) const; 
   void             FillLUT();
           Double_t LUTposition(Int_t ilayer, Double_t ampL, Double_t ampC, Double_t ampR) const;
-  virtual void     ResetHelperIndexes(AliTRDSignalIndex *indexesIn);
   
-  void              SetPadStatus(UChar_t status, UChar_t &encoding);
-  UChar_t           GetPadStatus(UChar_t encoding);
-  Int_t             GetCorruption(UChar_t encoding);
+  void             SetPadStatus(const UChar_t status, UChar_t &encoding);
+  UChar_t          GetPadStatus(UChar_t encoding) const;
+  Int_t            GetCorruption(UChar_t encoding) const;
+
+  Bool_t           IsMaximum(const Int_t row, const Int_t col, const Int_t time, 
+                            UChar_t &pasStatus, Double_t *const Signals);
+  Bool_t           IsFivePadCluster(const Int_t row, const Int_t col, const Int_t time, 
+                                   Double_t *SignalsThisMax, Double_t *SignalsNeighbourMax, Double_t &ratio);
+  void             CreateCluster(const Int_t row, const Int_t col, const Int_t time, 
+                                const Double_t *const clusterSignal, const UChar_t padStatus);
 
   const AliTRDReconstructor *fReconstructor;       //! reconstructor
   AliRunLoader        *fRunLoader;           //! Run Loader
@@ -115,14 +114,45 @@ class AliTRDclusterizer : public TNamed
   Bool_t               fAddLabels;           //  Should clusters have MC labels?
   Int_t                fRawVersion;          //  Expected raw version of the data - default is 2
 
-  AliTRDSignalIndex   *fIndexesOut;          //! Helper indexes for clusterization
-  AliTRDSignalIndex   *fIndexesMaxima;       //! Helper indexes for clusterization
-
   AliTRDtransform     *fTransform;           //! Transforms the reconstructed space points
 
   Int_t                fLUTbin;              //  Number of bins of the LUT
   Double_t            *fLUT;                 //! The lookup table
 
+  AliTRDarrayADC      *fDigitsIn;
+  AliTRDSignalIndex   *fIndexes;
+  Float_t              fADCthresh;            // ADC thresholds: There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO
+  Float_t              fMaxThresh;            // Threshold value for the maximum
+  Float_t              fSigThresh;            // Threshold value for the digit signal
+  Float_t              fMinMaxCutSigma;       // Threshold value for the maximum (cut noise)
+  Float_t              fMinLeftRightCutSigma; // Threshold value for the sum pad (cut noise)
+  Int_t                fLayer;                // Current layer of the detector
+  Int_t                fDet;                  // Current detecor
+  UShort_t             fVolid;                // Volume ID
+  Int_t                fColMax;               // Number of Colums in one detector
+  Int_t                fTimeTotal;            // Number of time bins
+  AliTRDCalROC        *fCalGainFactorROC;     // Calibration object with pad wise values for the gain factors
+  Float_t              fCalGainFactorDetValue;// Calibration value for chamber wise noise
+  AliTRDCalROC        *fCalNoiseROC;          // Calibration object with pad wise values for the noise
+  Float_t              fCalNoiseDetValue;     // Calibration value for chamber wise noise
+  AliTRDarraySignal   *fDigitsOut;
+  Int_t                fClusterROC;           // The index to the first cluster of a given ROC
+  Int_t                firstClusterROC;       // The number of cluster in a given ROC
+
+  struct ClusterizerStruct
+  {
+    Int_t       Row;
+    Int_t       Col;
+    Int_t       Time;
+    UChar_t     padStatus;
+    Double_t    Signals[3];
+    ClusterizerStruct():Row(0),Col(0),Time(0),padStatus(0)
+      {}
+    ClusterizerStruct &operator=(const ClusterizerStruct &a)
+      {Row=a.Row; Col=a.Col; Time=a.Time; padStatus=a.padStatus;
+       memcpy(Signals, a.Signals, 3*sizeof(Double_t)); return *this;}
+  };
+
   ClassDef(AliTRDclusterizer,6)              //  TRD clusterfinder
 
 };