Remove the last print statements
[u/mrichter/AliRoot.git] / TRD / AliTRDmcm.cxx
index 2e751da..481c81e 100644 (file)
 
 /*
 $Log$
+Revision 1.3  2006/08/11 17:58:05  cblume
+Next round of effc++ changes
+
+Revision 1.2  2006/04/05 12:45:40  hristov
+Updated TRD trigger code. Now the trigger code can run:
+
+- in the simulation step, through the central trigger interface, to
+produce and store "tracklets" and to produce and analyze tracks, with
+inputs to the central trigger
+
+- in the reconstruction step: if the tracklets have already been produced
+in the simulation step (by the trigger option), then only the tracks are
+produced and stored in the ESD event; if not, the trigger start by
+producing the tracklets, etc.
+
+Bogdan
+
 Revision 1.1.1.1  2004/08/19 14:58:11  vulpescu
 CVS head
 
@@ -23,8 +40,17 @@ test
 
 */
 
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//                                                                           //
+//  Multi Chip Module exponential filter and tracklet finder                 //
+//                                                                           //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
 #include <TMath.h>
-#include <TH2F.h>
+
+#include "AliLog.h"
 
 #include "AliTRDmcm.h"
 #include "AliTRDtrigParam.h"
@@ -33,99 +59,138 @@ ClassImp(AliTRDmcm)
 
 //_____________________________________________________________________________
 AliTRDmcm::AliTRDmcm() 
+  :TObject()
+  ,fTrigParam(0)
+  ,fNtrk(0)
+  ,fRobId(0)
+  ,fChaId(0)
+  ,fRow(0)
+  ,fColFirst(0)
+  ,fColLast(0)
+  ,fTime1(0)
+  ,fTime2(0)
+  ,fClusThr(0)
+  ,fPadThr(0)
+  ,fNtrkSeeds(0)
+  ,fR1(0)
+  ,fR2(0)
+  ,fC1(0)
+  ,fC2(0)
+  ,fPedestal(0)
+  ,fId(0)
 {
-
   //
   // AliTRDmcm default constructor
   //
 
-  fTrigParam = 0;
+  Int_t i = 0;
+  Int_t j = 0;
 
-  fNtrk = 0;
-  for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
+  for (i = 0; i < kMaxTrackletsPerMCM; i++) {
     fTrkIndex[i] = 0;
+    fSeedCol[i]  = -1;
   }
-  fRobId = 0;
-  fChaId = 0;
-  fRow = 0;
-  fColFirst = 0;
-  fColLast = 0;
-  for (Int_t i = 0; i < kMcmCol; i++) {
-    fPadHits[i] = 0;
-    for (Int_t j = 0; j < kMcmTBmax; j++) {
-      fADC[i][j] = 0.0;
+  for (i = 0; i < kMcmCol; i++) {
+    fPadHits[i]  = 0;
+    for (j = 0; j < kMcmTBmax; j++) {
+      fADC[i][j]    = 0.0;
       fIsClus[i][j] = kFALSE;
     }
   }
-  fPadThr = 0;
-  fClusThr = 0;
-  fTime1 = 0;
-  fTime2 = 0;
-  fNtrkSeeds = 0;
-  for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
-    fSeedCol[i] = -1;
-  }
-
-  fR1 = 0.0;
-  fR2 = 0.0;
-  fC1 = 0.0;
-  fC2 = 0.0;
-  fPedestal = 0.0;
-
-  fId = 0;
 
 }
 
 //_____________________________________________________________________________
 AliTRDmcm::AliTRDmcm(AliTRDtrigParam *trigp, const Int_t id) 
+  :TObject()
+  ,fTrigParam(trigp)
+  ,fNtrk(0)
+  ,fRobId(0)
+  ,fChaId(0)
+  ,fRow(0)
+  ,fColFirst(0)
+  ,fColLast(0)
+  ,fTime1(trigp->GetTime1())
+  ,fTime2(trigp->GetTime2())
+  ,fClusThr(trigp->GetClusThr())
+  ,fPadThr(trigp->GetPadThr())
+  ,fNtrkSeeds(0)
+  ,fR1(0)
+  ,fR2(0)
+  ,fC1(0)
+  ,fC2(0)
+  ,fPedestal(0)
+  ,fId(id)
 {
   //
   // AliTRDmcm constructor
   //
 
-  fTrigParam = trigp;
+  Int_t i = 0;
+  Int_t j = 0;
 
-  fNtrk = 0;
-  for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
+  for (i = 0; i < kMaxTrackletsPerMCM; i++) {
     fTrkIndex[i] = 0;
+    fSeedCol[i]  = -1;
   }
-  fRobId = 0;
-  fChaId = 0;
-  fRow = 0;
-  fColFirst = 0;
-  fColLast = 0;
-  for (Int_t i = 0; i < kMcmCol; i++) {
-    fPadHits[i] = 0;
-    for (Int_t j = 0; j < kMcmTBmax; j++) {
-      fADC[i][j] = 0.0;
+  for (i = 0; i < kMcmCol; i++) {
+    fPadHits[i]  = 0;
+    for (j = 0; j < kMcmTBmax; j++) {
+      fADC[i][j]    = 0.0;
       fIsClus[i][j] = kFALSE;
     }
   }
-  fPadThr  = fTrigParam->GetPadThr();
-  fClusThr = fTrigParam->GetClusThr();
-  fTime1 = fTrigParam->GetTime1();
-  fTime2 = fTrigParam->GetTime2();
-  fNtrkSeeds = 0;
-  for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
-    fSeedCol[i] = -1;
-  }
   
-  fR1 = 0.0;
-  fR2 = 0.0;
-  fC1 = 0.0;
-  fC2 = 0.0;
-  fPedestal = 0.0;
-
   fTrigParam->GetFilterParam(fR1,fR2,fC1,fC2,fPedestal);
 
-  fId = id;
+}
+
+//_____________________________________________________________________________
+AliTRDmcm::AliTRDmcm(const AliTRDmcm &m) 
+  :TObject(m)
+  ,fTrigParam(NULL)
+  ,fNtrk(m.fNtrk)
+  ,fRobId(m.fRobId)
+  ,fChaId(m.fChaId)
+  ,fRow(m.fRow)
+  ,fColFirst(m.fColFirst)
+  ,fColLast(m.fColLast)
+  ,fTime1(m.fTime1)
+  ,fTime2(m.fTime2)
+  ,fClusThr(m.fClusThr)
+  ,fPadThr(m.fPadThr)
+  ,fNtrkSeeds(m.fNtrkSeeds)
+  ,fR1(m.fR1)
+  ,fR2(m.fR2)
+  ,fC1(m.fC1)
+  ,fC2(m.fC2)
+  ,fPedestal(m.fPedestal)
+  ,fId(m.fId)
+{
+  //
+  // AliTRDmcm copy constructor
+  //
+
+  Int_t i = 0;
+  Int_t j = 0;
+
+  for (i = 0; i < kMaxTrackletsPerMCM; i++) {
+    ((AliTRDmcm &) m).fTrkIndex[i] = 0;
+    ((AliTRDmcm &) m).fSeedCol[i]  = -1;
+  }
+  for (i = 0; i < kMcmCol; i++) {
+    ((AliTRDmcm &) m).fPadHits[i]  = 0;
+    for (j = 0; j < kMcmTBmax; j++) {
+      ((AliTRDmcm &) m).fADC[i][j]    = 0.0;
+      ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
+    }
+  }
 
 }
 
 //_____________________________________________________________________________
 AliTRDmcm::~AliTRDmcm() 
 {
-
   //
   // AliTRDmcm destructor
   //
@@ -133,6 +198,62 @@ AliTRDmcm::~AliTRDmcm()
 }
 
 //_____________________________________________________________________________
+AliTRDmcm &AliTRDmcm::operator=(const AliTRDmcm &m)
+{
+  //
+  // Assignment operator
+  //
+
+  if (this != &m) ((AliTRDmcm &) m).Copy(*this); 
+
+  return *this;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDmcm::Copy(TObject &m) const
+{
+  //
+  // Copy function
+  //
+
+  Int_t i = 0;
+  Int_t j = 0;
+
+  ((AliTRDmcm &) m).fTrigParam  = NULL;
+  ((AliTRDmcm &) m).fNtrk       = fNtrk;
+  ((AliTRDmcm &) m).fRobId      = fRobId;
+  ((AliTRDmcm &) m).fChaId      = fChaId;
+  ((AliTRDmcm &) m).fRow        = fRow;
+  ((AliTRDmcm &) m).fColFirst   = fColFirst;
+  ((AliTRDmcm &) m).fColLast    = fColLast;
+  ((AliTRDmcm &) m).fTime1      = fTime1;
+  ((AliTRDmcm &) m).fTime2      = fTime2;
+  ((AliTRDmcm &) m).fClusThr    = fClusThr;
+  ((AliTRDmcm &) m).fPadThr     = fPadThr;
+  ((AliTRDmcm &) m).fNtrkSeeds  = fNtrkSeeds;
+  ((AliTRDmcm &) m).fR1         = fR1;
+  ((AliTRDmcm &) m).fR2         = fR2;
+  ((AliTRDmcm &) m).fC1         = fC1;
+  ((AliTRDmcm &) m).fC2         = fC2;
+  ((AliTRDmcm &) m).fPedestal   = fPedestal;
+  ((AliTRDmcm &) m).fId         = fId;
+
+  for (i = 0; i < kMaxTrackletsPerMCM; i++) {
+    ((AliTRDmcm &) m).fTrkIndex[i] = 0;
+    ((AliTRDmcm &) m).fSeedCol[i]  = -1;
+  }
+  for (i = 0; i < kMcmCol; i++) {
+    ((AliTRDmcm &) m).fPadHits[i]  = 0;
+    for (j = 0; j < kMcmTBmax; j++) {
+      ((AliTRDmcm &) m).fADC[i][j]    = 0.0;
+      ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
+    }
+  }
+
+}
+
+//_____________________________________________________________________________
 void AliTRDmcm::AddTrk(const Int_t id) 
 {
   //
@@ -153,14 +274,17 @@ void AliTRDmcm::Reset()
   // Reset MCM data
   //
 
-  for (Int_t i = 0; i < kMcmCol; i++) {
+  Int_t i = 0;
+  Int_t j = 0;
+
+  for (i = 0; i < kMcmCol; i++) {
     fPadHits[i] = 0;
-    for (Int_t j = 0; j < kMcmTBmax; j++) {
-      fADC[i][j] = 0.0;
+    for (j = 0; j < kMcmTBmax; j++) {
+      fADC[i][j]    = 0.0;
       fIsClus[i][j] = kFALSE;
     }
   }
-  for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
+  for (i = 0; i < kMaxTrackletsPerMCM; i++) {
     fSeedCol[i] = -1;
   }
   
@@ -173,87 +297,77 @@ Bool_t AliTRDmcm::Run()
   // Run MCM
   //
 
-  if ( fTrigParam->GetDebugLevel() > 1 ) printf("AliTRDmcm::Run MCM %d\n",Id());
+  AliDebug(2,(Form("Run MCM %d\n",Id())));
+
+  Int_t   iTime = 0;
+  Int_t   iCol  = 0;
+  Int_t   iPad  = 0;
+  Int_t   iPlus = 0;
+  Int_t   i     = 0;
+  Int_t   j     = 0;
 
-  Float_t Amp[3] = {0.0, 0.0, 0.0};
+  Float_t amp[3] = { 0.0, 0.0, 0.0 };
   Int_t   nClus;
-  Int_t   ClusCol[kMcmCol/2];
-  Float_t ClusAmp[kMcmCol/2];
-  Float_t VeryLarge;
-  Int_t   ClusMin = -1;
+  Int_t   clusCol[kMcmCol/2];
+  Float_t clusAmp[kMcmCol/2];
+  Float_t veryLarge;
+  Int_t   clusMin = -1;
   
-  for (Int_t iTime = fTime1; iTime <= fTime2; iTime++) {  // main TB loop
+  // Main TB loop
+  for (iTime = fTime1; iTime <= fTime2; iTime++) {  
 
-    // find clusters...
+    // Find clusters...
     nClus = 0;
-    for (Int_t iCol = 1; iCol < (kMcmCol-1); iCol++) {
-      Amp[0] = fADC[iCol-1][iTime];
-      Amp[1] = fADC[iCol  ][iTime];
-      Amp[2] = fADC[iCol+1][iTime];
-      if (IsCluster(Amp)) {
+    for (iCol = 1; iCol < (kMcmCol-1); iCol++) {
+      amp[0] = fADC[iCol-1][iTime];
+      amp[1] = fADC[iCol  ][iTime];
+      amp[2] = fADC[iCol+1][iTime];
+      if (IsCluster(amp)) {
        fIsClus[iCol][iTime] = kTRUE;
-       ClusCol[nClus] = iCol;
-       ClusAmp[nClus] = Amp[0]+Amp[1]+Amp[2];
+       clusCol[nClus]       = iCol;
+       clusAmp[nClus]       = amp[0]+amp[1]+amp[2];
        nClus++;
        if (nClus == kMcmCol/2) {
-         printf("Too many clusters in time bin %2d MCM %d...\n",iTime,Id());
-         //return kFALSE;
+         AliWarning(Form("Too many clusters in time bin %2d MCM %d...\n",iTime,Id()));
          break;
        }
       }
     }
 
     // ...but no more than six...
-    if (nClus > (Int_t)kSelClus) {
-      for (Int_t j = kSelClus/2; j < nClus-kSelClus/2; j++) {
-       fIsClus[ClusCol[j]][iTime] = kFALSE;
+    if (nClus > (Int_t) kSelClus) {
+      for (j = kSelClus/2; j < nClus-kSelClus/2; j++) {
+       fIsClus[clusCol[j]][iTime] = kFALSE;
       } 
     }
 
     // ...and take the largest four.
-
     Int_t nClusPlus = nClus - kMaxClus;
-    for (Int_t iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
-      VeryLarge = 1.E+10;
-      for (Int_t i = 0; i < nClus; i++) {
-       if (fIsClus[ClusCol[i]][iTime]) {
-         if (ClusAmp[i] <= VeryLarge) {
-           VeryLarge = ClusAmp[i];
-           ClusMin = i;
+    for (iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
+      veryLarge = 1.E+10;
+      for (i = 0; i < nClus; i++) {
+       if (fIsClus[clusCol[i]][iTime]) {
+         if (clusAmp[i] <= veryLarge) {
+           veryLarge = clusAmp[i];
+           clusMin = i;
          }
        }
       }
-      fIsClus[ClusCol[ClusMin]][iTime] = kFALSE;
+      fIsClus[clusCol[clusMin]][iTime] = kFALSE;
     }
 
     AddTimeBin(iTime);
 
   }  // end main TB loop
-  
-  if (fTrigParam->GetDebugLevel() > 1) {
-    for (Int_t i = fTime1; i <= fTime2; i++) {
-      printf("%2d: ",i);
-      for (Int_t j = 0; j < kMcmCol; j++) {
-       printf("%1d ",fIsClus[j][i]);
-      }
-      printf("\n");
-    }
-    printf("PadHits:     ");
-    for (Int_t iPad = 0; iPad < kMcmCol; iPad++) {
-      printf("%2d ",fPadHits[iPad]);
-    }
-    printf("\n");
-  }
-  
+    
   if ((fNtrkSeeds = CreateSeeds())) {
-
     return kTRUE;
-
   }
 
   return kFALSE;
 
 }
+
 //_____________________________________________________________________________
 Int_t AliTRDmcm::CreateSeeds()
 {
@@ -261,17 +375,19 @@ Int_t AliTRDmcm::CreateSeeds()
   // Make column seeds (from Falk Lesser, ex KIP)
   //
 
-  if ( fTrigParam->GetDebugLevel() > 1 ) printf("AliTRDmcm::CreateSeeds MCM %d \n",Id());
-
+  Int_t i      = 0;
+  Int_t j      = 0;
   Int_t nSeeds = 0;
 
-  // working array for hit sums
+  AliDebug(2,Form("AliTRDmcm::CreateSeeds MCM %d \n",Id()));
+
+  // Working array for hit sums
   Int_t fHit2padSum[2][kMcmCol];            
 
-  // initialize the array
-  for( Int_t i = 0; i < 2; i++ ) {
-    for( Int_t j = 0; j < kMcmCol; j++ ) {
-      if( i == 0 ) {
+  // Initialize the array
+  for (i = 0; i < 2; i++) {
+    for (j = 0; j < kMcmCol; j++ ) {
+      if (i == 0) {
        fHit2padSum[i][j] = j; 
       } else {
        fHit2padSum[i][j] = -1; 
@@ -279,106 +395,66 @@ Int_t AliTRDmcm::CreateSeeds()
     }
   }
 
-  Int_t Sum10 = fTrigParam->GetSum10();
-  Int_t Sum12 = fTrigParam->GetSum12();
+  Int_t sum10 = fTrigParam->GetSum10();
+  Int_t sum12 = fTrigParam->GetSum12();
 
-  // build the 2padSum
-  Int_t Nsum2seed = 0;
-  for( Int_t i = 0; i < kMcmCol; i++ ) {
-    if( i < (kMcmCol-1) ) {
-      if( (fPadHits[i] >= Sum10) && ((fPadHits[i] + fPadHits[i+1]) >= Sum12) ) {
+  // Build the 2padSum
+  Int_t nsum2seed = 0;
+  for (i = 0; i < kMcmCol; i++) {
+    if (i < (kMcmCol-1)) {
+      if ((fPadHits[i] >= sum10) && ((fPadHits[i] + fPadHits[i+1]) >= sum12)) {
        fHit2padSum[1][i] = fPadHits[i] + fPadHits[i+1]; 
-      } else {
+      } 
+      else {
        fHit2padSum[1][i] = -1;
       }
-    } else {
-      if ( fPadHits[i] >= Sum12 ) {
+    } 
+    else {
+      if (fPadHits[i] >= sum12) {
        fHit2padSum[1][i] = fPadHits[i]; 
-      } else {
+      } 
+      else {
        fHit2padSum[1][i] = -1;
       }
     }
-    if (fHit2padSum[1][i] > 0) Nsum2seed++;
-  }
-
-  if (fTrigParam->GetDebugLevel() > 1) {
-    printf("fHit2padSum: ");
-    for( Int_t i = 0; i < kMcmCol; i++ ) {
-      printf("%2d ",fHit2padSum[0][i]);
+    if (fHit2padSum[1][i] > 0) {
+      nsum2seed++;
     }
-    printf("\n");
-    printf("             ");
-    for( Int_t i = 0; i < kMcmCol; i++ ) {
-      printf("%2d ",fHit2padSum[1][i]);
-    }
-    printf("\n");
   }
 
   // sort the sums in decreasing order of the amplitude        
   Sort(kMcmCol,&fHit2padSum[0][0],&fHit2padSum[1][0],1);
 
-  if (fTrigParam->GetDebugLevel() > 1) {
-    printf("fHit2padSum: ");
-    for( Int_t i = 0; i < kMcmCol; i++ ) {
-      printf("%2d ",fHit2padSum[0][i]);
-    }
-    printf("\n");
-    printf("             ");
-    for( Int_t i = 0; i < kMcmCol; i++ ) {
-      printf("%2d ",fHit2padSum[1][i]);
-    }
-    printf("\n");
-  }
-
   // arrange (maximum number of) candidates in increasing order of the column number
-  nSeeds = TMath::Min(Nsum2seed,kMaxTrackletsPerMCM);
+  nSeeds = TMath::Min(nsum2seed,kMaxTrackletsPerMCM);
   Sort(nSeeds,&fHit2padSum[1][0],&fHit2padSum[0][0],0);
 
-  for (Int_t i = 0; i < nSeeds; i++) {
+  for (i = 0; i < nSeeds; i++) {
     fSeedCol[i] = fHit2padSum[0][i];
   }
 
-  if (fTrigParam->GetDebugLevel() > 1) {
-    printf("Found %d seeds before multiple rejection. \n",nSeeds);
-    printf("fHit2padSum: ");
-    for( Int_t i = 0; i < kMcmCol; i++ ) {
-      printf("%2d ",fHit2padSum[0][i]);
-    }
-    printf("\n");
-    printf("             ");
-    for( Int_t i = 0; i < kMcmCol; i++ ) {
-      printf("%2d ",fHit2padSum[1][i]);
-    }
-    printf("\n");
-  }
-
   // reject multiple found tracklets
-  Int_t imax = nSeeds-1;
-  for (Int_t i = 0; i < imax; i++) {
+  Int_t imax = nSeeds - 1;
+  for (i = 0; i < imax; i++) {
 
     if ((fHit2padSum[0][i]+1) == fHit2padSum[0][i+1]) {
       nSeeds--;
       if (fHit2padSum[1][i] >= fHit2padSum[1][i+1]) {
-       if (fTrigParam->GetDebugLevel() > 1) 
-         printf("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]);
+       if (fTrigParam->GetDebugLevel() > 1) {
+         AliWarning(Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]));
+       }
        fSeedCol[i+1] = -1;
-      } else {
-       if (fTrigParam->GetDebugLevel() > 1) 
-         printf("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]);
-       fSeedCol[i] = -1;
+      } 
+      else {
+       if (fTrigParam->GetDebugLevel() > 1) {
+         AliWarning(Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]));
+       }
+       fSeedCol[i]   = -1;
       }
     }
 
   }
 
-  if ( fTrigParam->GetDebugLevel() > 1 ) {
-    printf("Found %d seeds in MCM %d ",nSeeds,Id());
-    for (Int_t i = 0; i < (imax+1); i++) {
-      if (fSeedCol[i] >= 0) printf(", %02d ",fSeedCol[i]);
-    }
-    printf("\n");
-  }
-
   return nSeeds;
 
 }
@@ -386,47 +462,52 @@ Int_t AliTRDmcm::CreateSeeds()
 //_____________________________________________________________________________
 void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir)
 {
-
+  //
   // Sort two parallel vectors (x1[nel], x2[nel]) after the second one (x2)
   // in the direction: dir = 0 ascending order
   //                   dir = 1 descending order
+  //
 
+  Int_t  i = 0;
   Bool_t sort;
-  Int_t tmp1, tmp2;
+  Int_t  tmp1;
+  Int_t  tmp2;
 
-  if ( dir == 0 ) {
+  if (dir == 0) {
 
     do { 
       sort = kTRUE;
-      for ( Int_t i = 0; i < (nel-1); i++ )
-       if ( x2[i+1] < x2[i] ) {
-         tmp2 = x2[i]; 
-         x2[i] = x2[i+1]; 
+      for (i = 0; i < (nel-1); i++) {
+       if (x2[i+1] < x2[i]) {
+         tmp2    = x2[i]; 
+         x2[i]   = x2[i+1]; 
          x2[i+1] = tmp2;
-         tmp1 = x1[i]; 
-         x1[i] = x1[i+1]; 
+         tmp1    = x1[i]; 
+         x1[i]   = x1[i+1]; 
          x1[i+1] = tmp1;
-         sort = kFALSE;
+         sort    = kFALSE;
        }
-    } while ( sort == kFALSE );
+      }
+    } while (sort == kFALSE);
 
   }
 
-  if ( dir == 1 ) {
+  if (dir == 1) {
 
     do { 
       sort = kTRUE;
-      for ( Int_t i = 0; i < (nel-1); i++ )
-       if ( x2[i+1] > x2[i] ) {
-         tmp2 = x2[i]; 
-         x2[i] = x2[i+1]; 
+      for (i = 0; i < (nel-1); i++) {
+       if (x2[i+1] > x2[i]) {
+         tmp2    = x2[i]; 
+         x2[i]   = x2[i+1]; 
          x2[i+1] = tmp2;
-         tmp1 = x1[i]; 
-         x1[i] = x1[i+1]; 
+         tmp1    = x1[i]; 
+         x1[i]   = x1[i+1]; 
          x1[i+1] = tmp1;
-         sort = kFALSE;
+         sort    = kFALSE;
        }
-    } while ( sort == kFALSE );
+      }
+    } while (sort == kFALSE);
 
   }
 
@@ -448,20 +529,26 @@ void AliTRDmcm::AddTimeBin(const Int_t iTime)
 }
 
 //_____________________________________________________________________________
-Bool_t AliTRDmcm::IsCluster(Float_t amp[3])
+Bool_t AliTRDmcm::IsCluster(Float_t amp[3]) const
 {
   //
   // Find if the amplitudes amp[0], amp[1], amp[2] are a cluster
   //
 
   // -> shape
-  if (amp[0] > amp[1] || amp[2] > amp[1]) return kFALSE;
+  if (amp[0] > amp[1] || amp[2] > amp[1]) {
+    return kFALSE;
+  }
 
   // -> cluster amplitude
-  if ((amp[0]+amp[1]+amp[2]) < fClusThr) return kFALSE;
+  if ((amp[0]+amp[1]+amp[2]) < fClusThr) {
+    return kFALSE;
+  }
 
   // -> pad amplitude
-  if (amp[0] < fPadThr && amp[2] < fPadThr) return kFALSE;
+  if (amp[0] < fPadThr && amp[2] < fPadThr) {
+    return kFALSE;
+  }
 
   return kTRUE;
 
@@ -471,59 +558,60 @@ Bool_t AliTRDmcm::IsCluster(Float_t amp[3])
 void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
 {
   //
-  // exponential filter
+  // Exponential filter
   //
 
+  Int_t iCol  = 0;
+  Int_t iTime = 0;
+
   Double_t sour[kMcmTBmax];
-  Double_t Dtarg[kMcmTBmax];
-  Int_t    Itarg[kMcmTBmax];
+  Double_t dtarg[kMcmTBmax];
+  Int_t    itarg[kMcmTBmax];
 
   switch(ftype) {
 
-  case 0:
+   case 0:
 
-    for (Int_t iCol = 0; iCol < kMcmCol; iCol++) {
-      for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
+    for (iCol = 0; iCol < kMcmCol; iCol++) {
+      for (iTime = 0; iTime < kMcmTBmax; iTime++) {
        sour[iTime] = fADC[iCol][iTime];
       }
-      DeConvExpA(sour,Dtarg,kMcmTBmax,nexp);
-      for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
-       //fADC[iCol][iTime] = (Int_t)TMath::Max(0.0,Dtarg[iTime]);
-       fADC[iCol][iTime] = TMath::Max(0.0,Dtarg[iTime]);
+      DeConvExpA(sour,dtarg,kMcmTBmax,nexp);
+      for (iTime = 0; iTime < kMcmTBmax; iTime++) {
+       fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
       }
     }
     break;
 
   case 1:
 
-    for (Int_t iCol = 0; iCol < kMcmCol; iCol++) {
-      for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
+    for (iCol = 0; iCol < kMcmCol; iCol++) {
+      for (iTime = 0; iTime < kMcmTBmax; iTime++) {
        sour[iTime] = fADC[iCol][iTime];
       }
-      DeConvExpD(sour,Itarg,kMcmTBmax,nexp);
-      for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
-       fADC[iCol][iTime] = Itarg[iTime];
+      DeConvExpD(sour,itarg,kMcmTBmax,nexp);
+      for (iTime = 0; iTime < kMcmTBmax; iTime++) {
+       fADC[iCol][iTime] = itarg[iTime];
       }
     }
     break;
 
   case 2:
 
-    for (Int_t iCol = 0; iCol < kMcmCol; iCol++) {
-      for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
+    for (iCol = 0; iCol < kMcmCol; iCol++) {
+      for (iTime = 0; iTime < kMcmTBmax; iTime++) {
        sour[iTime] = fADC[iCol][iTime];
       }
-      DeConvExpMI(sour,Dtarg,kMcmTBmax);
-      for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
-       //fADC[iCol][iTime] = (Int_t)TMath::Max(0.0,Dtarg[iTime]);
-       fADC[iCol][iTime] = TMath::Max(0.0,Dtarg[iTime]);
+      DeConvExpMI(sour,dtarg,kMcmTBmax);
+      for (iTime = 0; iTime < kMcmTBmax; iTime++) {
+       fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
       }
     }
     break;
 
   default:
 
-    printf("Invalid filter type %d ! \n",ftype);
+    AliError(Form("Invalid filter type %d ! \n",ftype));
     return;
 
   }
@@ -533,270 +621,316 @@ void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
 //_____________________________________________________________________________
 void AliTRDmcm::DeConvExpA(Double_t *source, Double_t *target, Int_t n, Int_t nexp) 
 {
+  //
+  // Exponential filter "analog"
+  //
 
+  Int_t    i = 0;
+  Int_t    k = 0;
+  Double_t reminder[2];
+  Double_t correction;
+  Double_t result;
   Double_t rates[2];
   Double_t coefficients[2];
 
-  // initialize (coefficient = alpha, rates = lambda)
+  // Initialize (coefficient = alpha, rates = lambda)
   
   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
-  Double_t R1, R2, C1, C2;
-  R1 = (Double_t)fR1;
-  R2 = (Double_t)fR2;
-  C1 = (Double_t)fC1;
-  C2 = (Double_t)fC2;
+  Double_t r1, r2, c1, c2;
+  r1 = (Double_t)fR1;
+  r2 = (Double_t)fR2;
+  c1 = (Double_t)fC1;
+  c2 = (Double_t)fC2;
   
-  coefficients[0] = C1;
-  coefficients[1] = C2;
+  coefficients[0] = c1;
+  coefficients[1] = c2;
 
-  Double_t Dt = 0.100;
-  rates[0] = TMath::Exp(-Dt/(R1));
-  rates[1] = TMath::Exp(-Dt/(R2));
-  
-  Int_t i, k;
-  Double_t reminder[2];
-  Double_t correction, result;
+  Double_t dt = 0.1;
+  rates[0] = TMath::Exp(-dt/(r1));
+  rates[1] = TMath::Exp(-dt/(r2));
 
-  /* attention: computation order is important */
-  correction=0.0;
-  for ( k=0; k<nexp; k++ ) reminder[k]=0.0;
+  // Attention: computation order is important
+  correction = 0.0;
+  for (k = 0; k < nexp; k++) {
+    reminder[k] = 0.0;
+  }
     
-  for ( i=0; i<n; i++ ) {
-    result = ( source[i] - correction );    // no rescaling
+  for (i = 0; i < n; i++) {
+
+    result    = (source[i] - correction);    // no rescaling
     target[i] = result;
     
-    for ( k=0; k<nexp; k++ ) reminder[k] = rates[k] * ( reminder[k] + coefficients[k] * result);
+    for (k = 0; k < nexp; k++) {
+      reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
+    }
       
-    correction=0.0;
-    for ( k=0; k<nexp; k++ ) correction += reminder[k];
+    correction = 0.0;
+    for (k = 0; k < nexp; k++) {
+      correction += reminder[k];
+    }
+
   }
   
 }
 
 //_____________________________________________________________________________
-void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp) {
+void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp) 
+{
+  //
+  // Exponential filter "digital"
+  //
 
-  Int_t fAlpha_l, fAlpha_s, fLambda_l, fLambda_s, fTailPed;
-  Int_t iAlpha_l, iAlpha_s, iLambda_l, iLambda_s;
+  Int_t i = 0;
 
-  // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
-  // initialize (coefficient = alpha, rates = lambda)
-
-  fLambda_l = 0; fAlpha_l = 0; fLambda_s = 0; fAlpha_s = 0;
-  iLambda_l = 0; iAlpha_l = 0; iLambda_s = 0; iAlpha_s = 0;
+  Int_t fAlphaL;
+  Int_t fAlphaS;
+  Int_t fLambdaL;
+  Int_t fLambdaS;
+  Int_t fTailPed;
 
-  Double_t Dt = 0.100;
+  Int_t iAlphaL;
+  Int_t iAlphaS;
+  Int_t iLambdaL;
+  Int_t iLambdaS;
 
-  Double_t R1, R2, C1, C2;
-  R1 = (Double_t)fR1;
-  R2 = (Double_t)fR2;
-  C1 = (Double_t)fC1;
-  C2 = (Double_t)fC2;
+  // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
+  // initialize (coefficient = alpha, rates = lambda)
 
-  fLambda_l = (Int_t)((TMath::Exp(-Dt/R1)-0.75)*2048.0);
-  fLambda_s = (Int_t)((TMath::Exp(-Dt/R2)-0.25)*2048.0);
-  iLambda_l = fLambda_l & 0x01FF; iLambda_l |= 0x0600;         //  9 bit paramter + fixed bits
-  iLambda_s = fLambda_s & 0x01FF; iLambda_s |= 0x0200;         //  9 bit paramter + fixed bits
+  fLambdaL = 0; 
+  fAlphaL  = 0; 
+  fLambdaS = 0; 
+  fAlphaS  = 0;
+  iLambdaL = 0; 
+  iAlphaL  = 0; 
+  iLambdaS = 0; 
+  iAlphaS  = 0;
+
+  Double_t dt = 0.1;
+
+  Double_t r1;
+  Double_t r2;
+  Double_t c1;
+  Double_t c2;
+  r1 = (Double_t) fR1;
+  r2 = (Double_t) fR2;
+  c1 = (Double_t) fC1;
+  c2 = (Double_t) fC2;
+
+  fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
+  fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
+  iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600;    //  9 bit paramter + fixed bits
+  iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200;    //  9 bit paramter + fixed bits
 
   if (nexp == 1) {
-    fAlpha_l = (Int_t)(C1*2048.0);
-    iAlpha_l = fAlpha_l & 0x03FF;                              // 10 bit paramter
+    fAlphaL = (Int_t) (c1 * 2048.0);
+    iAlphaL = fAlphaL & 0x03FF;                                // 10 bit paramter
   }
   if (nexp == 2) {
-    fAlpha_l = (Int_t)(C1*2048.0);
-    fAlpha_s = (Int_t)((C2-0.5)*2048.0);
-    iAlpha_l = fAlpha_l & 0x03FF;                              // 10 bit paramter
-    iAlpha_s = fAlpha_s & 0x03FF; iAlpha_s |= 0x0400;          // 10 bit paramter + fixed bits
+    fAlphaL = (Int_t) (c1 * 2048.0);
+    fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
+    iAlphaL = fAlphaL & 0x03FF;                                // 10 bit paramter
+    iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400;             // 10 bit paramter + fixed bits
   }
   
-  Double_t iAl = iAlpha_l / 2048.0;           // alpha L: correspondence to floating point numbers
-  Double_t iAs = iAlpha_s / 2048.0;           // alpha S: correspondence to floating point numbers
-  Double_t iLl = iLambda_l / 2048.0;          // lambda L: correspondence to floating point numbers
-  Double_t iLs = iLambda_s / 2048.0;          // lambda S: correspondence to floating point numbers
-
-  Int_t h1,h2;
-  Int_t rem1, rem2;
-  Int_t correction, result;
-  Int_t iFactor = ((Int_t)fPedestal)<<2;
-
-  Double_t xi = 1-(iLl*iAs + iLs*iAl);             // calculation of equilibrium values of the
-  rem1 = (Int_t)((iFactor/xi) * ((1-iLs)*iLl*iAl)); // internal registers to prevent switch on effects.
-  rem2 = (Int_t)((iFactor/xi) * ((1-iLl)*iLs*iAs));
+  Double_t iAl = iAlphaL  / 2048.0;           // alpha L: correspondence to floating point numbers
+  Double_t iAs = iAlphaS  / 2048.0;           // alpha S: correspondence to floating point numbers
+  Double_t iLl = iLambdaL / 2048.0;           // lambda L: correspondence to floating point numbers
+  Double_t iLs = iLambdaS / 2048.0;           // lambda S: correspondence to floating point numbers
+
+  Int_t h1;
+  Int_t h2;
+  Int_t rem1;
+  Int_t rem2;
+  Int_t correction;
+  Int_t result;
+  Int_t iFactor = ((Int_t) fPedestal) << 2;
+
+  Double_t xi = 1 - (iLl*iAs + iLs*iAl);            // Calculation of equilibrium values of the
+  rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
+  rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
   
   // further initialization
-  if ( (rem1 + rem2) > 0x0FFF ) {
+  if ((rem1 + rem2) > 0x0FFF) {
     correction = 0x0FFF;
-  } else {
+  } 
+  else {
     correction = (rem1 + rem2) & 0x0FFF;
   }
 
   fTailPed = iFactor - correction;
 
-  for (Int_t i=0; i<n; i++) {
+  for (i = 0; i < n; i++) {
 
-    result = ( (Int_t)source[i] - correction );
-    if ( result<0 ) {                          
+    result = ((Int_t)source[i] - correction);
+    if (result < 0) {                          
       result = 0;
     }
 
     target[i] = result;
                                                         
-    h1 = ( rem1 + ((iAlpha_l * result) >> 11) );
-    if ( h1 > 0x0FFF ) {
+    h1 = (rem1 + ((iAlphaL * result) >> 11));
+    if (h1 > 0x0FFF) {
       h1 = 0x0FFF;
-    } else {
+    } 
+    else {
       h1 &= 0x0FFF;
     }
 
-    h2 = ( rem2 + ((iAlpha_s * result) >> 11));
-    if ( h2 > 0x0FFF ) {
+    h2 = (rem2 + ((iAlphaS * result) >> 11));
+    if (h2 > 0x0FFF) {
       h2 = 0x0FFF;
-    } else {
+    } 
+    else {
       h2 &= 0x0FFF;
     }
   
-    rem1 = (iLambda_l * h1 ) >> 11;
-    rem2 = (iLambda_s * h2 ) >> 11;
+    rem1 = (iLambdaL * h1 ) >> 11;
+    rem2 = (iLambdaS * h2 ) >> 11;
     
-    if ( (rem1 + rem2) > 0x0FFF ) {
+    if ((rem1 + rem2) > 0x0FFF) {
       correction = 0x0FFF;
-    } else {
+    } 
+    else {
       correction = (rem1 + rem2) & 0x0FFF;
     }
+
   }
 
 }
 
 //_____________________________________________________________________________
-void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n) {
+void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n) 
+{
+  //
+  // Exponential filter (M. Ivanov)
+  //
+
+  Int_t i = 0;
 
-   Double_t Sig1[100], Sig2[100], Sig3[100];//, Sig4[100];
-   for (Int_t i = 0; i < n; i++) Sig1[i] = source[i];
+  Double_t sig1[100];
+  Double_t sig2[100];
+  Double_t sig3[100];
 
-   Float_t Dt = 0.100;
+  for (i = 0; i < n; i++) {
+    sig1[i] = source[i];
+  }
 
-   //Float_t lambda0 = 9.8016*Dt;  // short
-   //Float_t lambda1 = 1.0778*Dt;  // long
+  Float_t dt = 0.1;
 
-   Float_t lambda0 = (1.0/fR2)*Dt;
-   Float_t lambda1 = (1.0/fR1)*Dt;
+  Float_t lambda0 = (1.0 / fR2) * dt;
+  Float_t lambda1 = (1.0 / fR1) * dt;
 
-   TailMakerSpline(Sig1,Sig2,lambda0,n);
-   TailCancelationMI(Sig2,Sig3,0.7,lambda1,n);
+  TailMakerSpline(sig1,sig2,lambda0,n);
+  TailCancelationMI(sig2,sig3,0.7,lambda1,n);
 
-   for (Int_t i = 0; i < n; i++) target[i] = Sig3[i];
+  for (i = 0; i < n; i++) {
+    target[i] = sig3[i];
+  }
 
 }
 
 //______________________________________________________________________________
-void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n) {
-
-   Double_t l = TMath::Exp(-lambda*0.5);
-   //
-   //
-   Double_t in[1000];
-   Double_t out[1000];
-   // initialize in[] and out[] goes 0 ... 2*n+19
-   for (Int_t i=0; i<n*2+20; i++){
-     in[i] = 0;
-     out[i]= 0;
-   }
-
-   // in[] goes 0, 1
-   in[0] = ampin[0];
-   in[1] = (ampin[0]+ampin[1])*0.5;
+void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n) 
+{
+  //
+  // Special filter (M. Ivanov)
+  //
+
+  Int_t    i = 0;
+
+  Double_t l = TMath::Exp(-lambda*0.5);
+  Double_t in[1000];
+  Double_t out[1000];
+
+  // Initialize in[] and out[] goes 0 ... 2*n+19
+  for (i = 0; i < n*2+20; i++) {
+    in[i]  = 0;
+    out[i] = 0;
+  }
+
+  // in[] goes 0, 1
+  in[0] = ampin[0];
+  in[1] = (ampin[0] + ampin[1]) * 0.5;
    
-   // add charge to the end
-   for (Int_t i=0; i<22; i++){
-     // in[] goes 2*n-2, 2*n-1, ... , 2*n+19 
-     in[2*(n-1)+i] = ampin[n-1];
-   }
-
-   // use arithm mean
-   for (Int_t i=1; i<n-1; i++){
-     // in[] goes 2, 3, ... , 2*n-4, 2*n-3
-     in[2*i]    = ampin[i];
-     in[2*i+1]  = ((ampin[i]+ampin[i+1]))/2.;
-   }
-   /*
-   // add charge to the end
-   for (Int_t i=-2; i<22; i++){
-     // in[] goes 2*n-4, 2*n-3, ... , 2*n+19 
-     in[2*(n-1)+i] = ampin[n-1];
-   }
-
-   // use spline mean
-   for (Int_t i=1; i<n-2; i++){
-     // in[] goes 2, 3, ... , 2*n-6, 2*n-5
-     in[2*i]    = ampin[i];
-     in[2*i+1]  = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16;
-   }
-   */
-   //
-   Double_t temp;
-   out[2*n]      = in[2*n];
-   temp          = 0;
-   for (int i=2*n; i>=0; i--){
-     out[i]      = in[i]   + temp;
-     temp        = l*(temp+in[i]);
-   }
-
-   //
-   for (int i=0;i<n;i++){
-     //ampout[i]      = out[2*i+1];  // org
-     ampout[i]      = out[2*i];
-   }
+  // Add charge to the end
+  for (i = 0; i < 22; i++) {
+    // in[] goes 2*n-2, 2*n-1, ... , 2*n+19 
+    in[2*(n-1)+i] = ampin[n-1];
+  }
+
+  // Use arithmetic mean
+  for (i = 1; i < n-1; i++) {
+    // in[] goes 2, 3, ... , 2*n-4, 2*n-3
+    in[2*i]   = ampin[i];
+    in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
+  }
+
+  Double_t temp;
+  out[2*n]    = in[2*n];
+  temp        = 0;
+  for (i = 2*n; i >= 0; i--) {
+    out[i]    = in[i] + temp;
+    temp      = l*(temp+in[i]);
+  }
+
+  for (i = 0; i < n; i++){
+    //ampout[i] = out[2*i+1];  // org
+    ampout[i] = out[2*i];
+  }
 
 }
 
 //______________________________________________________________________________
-void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout, Double_t norm, Double_t lambda, Int_t n) {
-
-   Double_t L = TMath::Exp(-lambda*0.5);
-   Double_t K = L*(1.-norm*lambda*0.5);
-   //
-   //
-   Double_t in[1000];
-   Double_t out[1000];
-   // initialize in[] and out[] goes 0 ... 2*n+19
-   for (Int_t i=0; i<n*2+20; i++){
-     in[i] = 0;
-     out[i]= 0;
-   }
-
-   // in[] goes 0, 1
-   in[0] = ampin[0];
-   in[1] = (ampin[0]+ampin[1])*0.5;
-
-   // add charge to the end
-   for (Int_t i=-2; i<22; i++){
-     // in[] goes 2*n-4, 2*n-3, ... , 2*n+19 
-     in[2*(n-1)+i] = ampin[n-1];
-   }
-
-   //
-   for (Int_t i=1; i<n-2; i++){
-     // in[] goes 2, 3, ... , 2*n-6, 2*n-5
-     in[2*i]    = ampin[i];
-     in[2*i+1]  = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16.;
-     //in[2*i+1]  = ((ampin[i]+ampin[i+1]))/2.;
-   }
-
-   //
-   Double_t temp;
-   out[0]     = in[0];
-   temp       = in[0];
-   for (int i=1; i<=2*n; i++){
-     out[i]      = in[i]   + (K-L)*temp;
-     temp        = in[i]   +  K*temp;
-   }
-   //
-   //
-   for (int i=0; i<n; i++){
-     //ampout[i]      = out[2*i+1];  // org
-     //ampout[i]      = TMath::Max(out[2*i+1], 0.0);  // org
-     ampout[i]      = TMath::Max(out[2*i], 0.0);
-   }
+void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout
+                                , Double_t norm, Double_t lambda, Int_t n) 
+{
+  //
+  // Special filter (M. Ivanov)
+  //
+
+  Int_t    i = 0;
+
+  Double_t l = TMath::Exp(-lambda*0.5);
+  Double_t k = l*(1.0 - norm*lambda*0.5);
+  Double_t in[1000];
+  Double_t out[1000];
+
+  // Initialize in[] and out[] goes 0 ... 2*n+19
+  for (i = 0; i < n*2+20; i++) {
+    in[i]  = 0;
+    out[i] = 0;
+  }
+
+  // in[] goes 0, 1
+  in[0] = ampin[0];
+  in[1] = (ampin[0]+ampin[1])*0.5;
+
+  // Add charge to the end
+  for (i =-2; i < 22; i++) {
+    // in[] goes 2*n-4, 2*n-3, ... , 2*n+19 
+    in[2*(n-1)+i] = ampin[n-1];
+  }
+
+  for (i = 1; i < n-2; i++) {
+    // in[] goes 2, 3, ... , 2*n-6, 2*n-5
+    in[2*i]    = ampin[i];
+    in[2*i+1]  = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
+    //in[2*i+1]  = ((ampin[i]+ampin[i+1]))/2.0;
+  }
+
+  Double_t temp;
+  out[0] = in[0];
+  temp   = in[0];
+  for (i = 1; i <= 2*n; i++) {
+    out[i] = in[i] + (k-l)*temp;
+    temp   = in[i] +  k   *temp;
+  }
+
+  for (i = 0; i < n; i++) {
+    //ampout[i] = out[2*i+1];  // org
+    //ampout[i] = TMath::Max(out[2*i+1],0.0);  // org
+    ampout[i] = TMath::Max(out[2*i],0.0);
+  }
 
 }