Fixes for #86059: Install data when ALICE_ROOT does not point to source (Christian)
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderMLEM.cxx
index 4fb83ca..21d26f5 100644 (file)
 #include <TMinuit.h>
 #include <TCanvas.h>
 #include <TMath.h>
-#include <TROOT.h>
 #include "AliCodeTimer.h"
 
 /// \cond CLASSIMP
 ClassImp(AliMUONClusterFinderMLEM)
 /// \endcond
  
-const Double_t AliMUONClusterFinderMLEM::fgkZeroSuppression = 6; // average zero suppression value
-//const Double_t AliMUONClusterFinderMLEM::fgkDistancePrecision = 1e-6; // (cm) used to check overlaps and so on
 const Double_t AliMUONClusterFinderMLEM::fgkDistancePrecision = 1e-3; // (cm) used to check overlaps and so on
 const TVector2 AliMUONClusterFinderMLEM::fgkIncreaseSize(-AliMUONClusterFinderMLEM::fgkDistancePrecision,-AliMUONClusterFinderMLEM::fgkDistancePrecision);
 const TVector2 AliMUONClusterFinderMLEM::fgkDecreaseSize(AliMUONClusterFinderMLEM::fgkDistancePrecision,AliMUONClusterFinderMLEM::fgkDistancePrecision);
 
- TMinuit* AliMUONClusterFinderMLEM::fgMinuit = 0x0;
 // Status flags for pads
 const Int_t AliMUONClusterFinderMLEM::fgkZero = 0x0; ///< pad "basic" state
 const Int_t AliMUONClusterFinderMLEM::fgkMustKeep = 0x1; ///< do not kill (for pixels)
@@ -76,21 +72,21 @@ fClusterList(),
 fEventNumber(0),
 fDetElemId(-1),
 fClusterNumber(0),
-fCathBeg(0),
+fHistMlem(0x0),
+fHistAnode(0x0),
 fPixArray(new TObjArray(20)),
 fDebug(0),
 fPlot(plot),
 fSplitter(0x0),
 fNClusters(0),
-fNAddVirtualPads(0)
+fNAddVirtualPads(0),
+fLowestPixelCharge(0),
+fLowestPadCharge(0),
+fLowestClusterCharge(0)
 {
   /// Constructor
-  fSegmentation[1] = fSegmentation[0] = 0x0; 
-
-  fPadBeg[0] = fPadBeg[1] = fCathBeg;
 
-  if (!fgMinuit) fgMinuit = new TMinuit(8);
+  fkSegmentation[1] = fkSegmentation[0] = 0x0; 
 
   if (fPlot) fDebug = 1;
 }
@@ -99,7 +95,7 @@ fNAddVirtualPads(0)
 AliMUONClusterFinderMLEM::~AliMUONClusterFinderMLEM()
 {
 /// Destructor
-  delete fgMinuit; fgMinuit = 0; delete fPixArray; fPixArray = 0;
+  delete fPixArray; fPixArray = 0;
 //  delete fDraw;
   delete fPreClusterFinder;
   delete fSplitter;
@@ -110,30 +106,36 @@ AliMUONClusterFinderMLEM::~AliMUONClusterFinderMLEM()
 //_____________________________________________________________________________
 Bool_t 
 AliMUONClusterFinderMLEM::Prepare(Int_t detElemId,
-                                  TClonesArray* pads[2],
+                                  TObjArray* pads[2],
                                   const AliMpArea& area,
                                   const AliMpVSegmentation* seg[2])
 {
   /// Prepare for clustering
-//  AliCodeTimerAuto("")
+//  AliCodeTimerAuto("",0)
   
   for ( Int_t i = 0; i < 2; ++i )
   {
-    fSegmentation[i] = seg[i];
+    fkSegmentation[i] = seg[i];
   }
   
   // Find out the DetElemId
   fDetElemId = detElemId;
   
   delete fSplitter;
-  fSplitter = new AliMUONClusterSplitterMLEM(fDetElemId,fPixArray);
+  fSplitter = new AliMUONClusterSplitterMLEM(fDetElemId,
+                                             fPixArray,
+                                             fLowestPixelCharge,
+                                             fLowestPadCharge,
+                                             fLowestClusterCharge);
   fSplitter->SetDebug(fDebug);
     
   // find out current event number, and reset the cluster number
-  fEventNumber = AliRunLoader::GetRunLoader()->GetEventNumber();
+  AliRunLoader *runLoader = AliRunLoader::Instance();
+  fEventNumber = runLoader ? runLoader->GetEventNumber() : 0;
   fClusterNumber = -1;
   fClusterList.Delete();
-  
+  fPixArray->Delete();
+
   AliDebug(3,Form("EVT %d DE %d",fEventNumber,fDetElemId));
   
   if ( fPreClusterFinder->NeedSegmentation() )
@@ -151,10 +153,17 @@ AliMUONCluster*
 AliMUONClusterFinderMLEM::NextCluster()
 {
   /// Return next cluster
-//  AliCodeTimerAuto("")
+//  AliCodeTimerAuto("",0)
   
   // if the list of clusters is not void, pick one from there
-  TObject* o = fClusterList.At(++fClusterNumber);
+  TObject* o(0x0);
+  
+  // do we have clusters in our list ?
+  if ( fClusterNumber < fClusterList.GetLast() )
+  {
+    o = fClusterList.At(++fClusterNumber);
+  }
+  
   if ( o != 0x0 ) return static_cast<AliMUONCluster*>(o);
   
   //FIXME : at this point, must check whether we've used all the digits
@@ -165,16 +174,17 @@ AliMUONClusterFinderMLEM::NextCluster()
   // pre-cluster and treat it
 
   fPreCluster = fPreClusterFinder->NextCluster();
-  
+
+  fPixArray->Delete();
+  fClusterList.Delete(); // reset the list of clusters for this pre-cluster
+  fClusterNumber = -1; //AZ
+    
   if (!fPreCluster)
   {
     // we are done
     return 0x0;
   }
     
-  fClusterList.Delete(); // reset the list of clusters for this pre-cluster
-  fClusterNumber = -1; //AZ
-  
   WorkOnPreCluster();
 
   // WorkOnPreCluster may have used only part of the pads, so we check that
@@ -200,10 +210,10 @@ AliMUONClusterFinderMLEM::WorkOnPreCluster()
   /// Starting from a precluster, builds a pixel array, and then
   /// extract clusters from this array
   
-  //  AliCodeTimerAuto("")     
+  //  AliCodeTimerAuto("",0)   
 
   if (fDebug) {
-    cout << " *** Event # " << AliRunLoader::GetRunLoader()->GetEventNumber() 
+    cout << " *** Event # " << fEventNumber 
         << " det. elem.: " << fDetElemId << endl;
     for (Int_t j = 0; j < fPreCluster->Multiplicity(); ++j) {
       AliMUONPad* pad = fPreCluster->Pad(j);
@@ -269,10 +279,9 @@ AliMUONClusterFinderMLEM::WorkOnPreCluster()
       }
     }
   } // for (Int_t i=0; i<nMax;
-  //if (nMax > 1) ((TH2D*) gROOT->FindObject("anode"))->Delete();
-  delete ((TH2D*) gROOT->FindObject("anode"));
-  //TH2D *mlem = (TH2D*) gROOT->FindObject("mlem");
-  //if (mlem) mlem->Delete();
+  delete fHistMlem;
+  delete fHistAnode;
+  fHistMlem = fHistAnode = 0x0;
   delete cluster;
   return kTRUE;
 }
@@ -298,11 +307,11 @@ AliMUONClusterFinderMLEM::CheckPrecluster(const AliMUONCluster& origCluster)
   /// Check precluster in order to attempt to simplify it (mostly for
   /// two-cathode preclusters)
     
-  AliCodeTimerAuto("")
+  AliCodeTimerAuto("",0)
 
   // Disregard small clusters (leftovers from splitting or noise)
   if ((origCluster.Multiplicity()==1 || origCluster.Multiplicity()==2) &&
-      origCluster.Charge(0)+origCluster.Charge(1) < 10) 
+      origCluster.Charge(0)+origCluster.Charge(1) < fLowestClusterCharge )
   { 
     return 0x0;
   }
@@ -368,10 +377,11 @@ AliMUONClusterFinderMLEM::CheckPreclusterTwoCathodes(AliMUONCluster* cluster)
       Int_t cath = pad->Cathode();
       Int_t cath1 = TMath::Even(cath);
       // Check for edge effect (missing pads on the _other_ cathode)
-      AliMpPad mpPad = fSegmentation[cath1]->PadByPosition(pad->Position(),kFALSE);
+      AliMpPad mpPad =
+      fkSegmentation[cath1]->PadByPosition(pad->Position().X(),
+                                           pad->Position().Y(),kFALSE);
       if (!mpPad.IsValid()) continue;
-      //if (nFlags == 1 && pad->Charge() < fgkZeroSuppression * 3) continue;
-      if (nFlags == 1 && pad->Charge() < 20) continue;
+      if (nFlags == 1 && pad->Charge() < fLowestPadCharge) continue; 
       AliDebug(2,Form("Releasing the following pad : de,cath,ix,iy %d,%d,%d,%d charge %e",
                       fDetElemId,pad->Cathode(),pad->Ix(),pad->Iy(),pad->Charge()));
       toBeRemoved.AddLast(pad);
@@ -541,7 +551,7 @@ AliMUONClusterFinderMLEM::CheckOverlaps()
                         cout << " Surface = " << pixelI->Size(0)*pixelI->Size(1)*4 << endl;
                         pixelJ->Print();
                         cout << " Surface = " << pixelJ->Size(0)*pixelJ->Size(1)*4 << endl;
-                        cout << " Area surface = " << area.Dimensions().X()*area.Dimensions().Y()*4 << endl;
+                        cout << " Area surface = " << area.GetDimensionX()*area.GetDimensionY()*4 << endl;
                         cout << "-------" << endl;
                         );
        */        
@@ -594,8 +604,8 @@ void AliMUONClusterFinderMLEM::BuildPixArrayOneCathode(AliMUONCluster& cluster)
 //  AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
 
   TVector2 dim = cluster.MinPadDimensions (-1, kFALSE);
-  Double_t width[2] = {dim.X(), dim.Y()}, xy0[2];
-  Int_t found[2] = {0}, mult = cluster.Multiplicity();
+  Double_t width[2] = {dim.X(), dim.Y()}, xy0[2]={99999,99999};
+  Int_t found[2] = {0,0}, mult = cluster.Multiplicity();
 
   for ( Int_t i = 0; i < mult; ++i) {
     AliMUONPad* pad = cluster.Pad(i);
@@ -613,24 +623,27 @@ void AliMUONClusterFinderMLEM::BuildPixArrayOneCathode(AliMUONCluster& cluster)
   if (cluster.Multiplicity(0) == 0) cath0 = 1;
   else if (cluster.Multiplicity(1) == 0) cath1 = 0;
 
-  TVector2 leftDown = cluster.Area(cath0).LeftDownCorner();
-  TVector2 rightUp = cluster.Area(cath0).RightUpCorner();
-  min[0] = leftDown.X();
-  min[1] = leftDown.Y();
-  max[0] = rightUp.X();
-  max[1] = rightUp.Y();
+
+  Double_t leftDownX, leftDownY;
+  cluster.Area(cath0).LeftDownCorner(leftDownX, leftDownY);
+  Double_t rightUpX, rightUpY;
+  cluster.Area(cath0).RightUpCorner(rightUpX, rightUpY);
+  min[0] = leftDownX;
+  min[1] = leftDownY;
+  max[0] = rightUpX;
+  max[1] = rightUpY;;
   if (cath1 != cath0) {
-    leftDown = cluster.Area(cath1).LeftDownCorner();
-    rightUp = cluster.Area(cath1).RightUpCorner();
-    min[0] = TMath::Max (min[0], leftDown.X());
-    min[1] = TMath::Max (min[1], leftDown.Y());
-    max[0] = TMath::Min (max[0], rightUp.X());
-    max[1] = TMath::Min (max[1], rightUp.Y());
+    cluster.Area(cath1).LeftDownCorner(leftDownX, leftDownY);
+    cluster.Area(cath1).RightUpCorner(rightUpX, rightUpY);
+    min[0] = TMath::Max (min[0], leftDownX);
+    min[1] = TMath::Max (min[1], leftDownY);
+    max[0] = TMath::Min (max[0], rightUpX);
+    max[1] = TMath::Min (max[1], rightUpY);
   }
 
   // Adjust limits
   //width[0] /= 2; width[1] /= 2; // just for check
-  Int_t nbins[2];
+  Int_t nbins[2]={0,0};
   for (Int_t i = 0; i < 2; ++i) {
     Double_t dist = (min[i] - xy0[i]) / width[i] / 2;
     if (TMath::Abs(dist) < 1.e-6) dist = -1.e-6;
@@ -653,7 +666,7 @@ void AliMUONClusterFinderMLEM::BuildPixArrayOneCathode(AliMUONCluster& cluster)
     AliMUONPad* pad = cluster.Pad(i);
     Int_t ix0 = xaxis->FindBin(pad->X());
     Int_t iy0 = yaxis->FindBin(pad->Y());
-    PadOverHist(0, ix0, iy0, pad);
+    PadOverHist(0, ix0, iy0, pad, hist1, hist2);
   }
 
   // Store pixels
@@ -691,12 +704,11 @@ void AliMUONClusterFinderMLEM::BuildPixArrayOneCathode(AliMUONCluster& cluster)
 }
 
 //_____________________________________________________________________________
-void AliMUONClusterFinderMLEM::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, AliMUONPad *pad)
+void AliMUONClusterFinderMLEM::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, AliMUONPad *pad,
+                                          TH2D *hist1, TH2D *hist2)
 {
   /// "Span" pad over histogram in the direction idir
 
-  TH2D *hist1 = static_cast<TH2D*> (gROOT->FindObject("Grid"));
-  TH2D *hist2 = static_cast<TH2D*> (gROOT->FindObject("Entries"));
   TAxis *axis = idir == 0 ? hist1->GetXaxis() : hist1->GetYaxis();
   Int_t nbins = axis->GetNbins(), cath = pad->Cathode();
   Double_t bin = axis->GetBinWidth(1), amask = TMath::Power(1000.,cath*1.);
@@ -708,7 +720,7 @@ void AliMUONClusterFinderMLEM::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, Ali
     if (ixy > nbins) break;
     Double_t lowEdge = axis->GetBinLowEdge(ixy);
     if (lowEdge + fgkDistancePrecision > pad->Coord(idir) + pad->Size(idir)) break;
-    if (idir == 0) PadOverHist(1, ixy, iy0, pad); // span in the other direction
+    if (idir == 0) PadOverHist(1, ixy, iy0, pad, hist1, hist2); // span in the other direction
     else {
       // Fill histogram
       Double_t cont = pad->Charge();
@@ -725,7 +737,7 @@ void AliMUONClusterFinderMLEM::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, Ali
     if (ixy < 1) break;
     Double_t upEdge = axis->GetBinUpEdge(ixy);
     if (upEdge - fgkDistancePrecision < pad->Coord(idir) - pad->Size(idir)) break;
-    if (idir == 0) PadOverHist(1, ixy, iy0, pad); // span in the other direction
+    if (idir == 0) PadOverHist(1, ixy, iy0, pad, hist1, hist2); // span in the other direction
     else {
       // Fill histogram
       Double_t cont = pad->Charge();
@@ -740,20 +752,20 @@ void AliMUONClusterFinderMLEM::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, Ali
 
 //_____________________________________________________________________________
 void
-AliMUONClusterFinderMLEM::Plot(const char* basename)
+AliMUONClusterFinderMLEM::Plot(const char* /*basename*/)
 {
   /// Make a plot and save it as png
   
   return; //AZ
-  if (!fPlot) return;
-  
-  TCanvas* c = new TCanvas("MLEM","MLEM",800,600);
-  c->Draw();
-  Draw();
-  c->Modified();
-  c->Update();
-  c->Print(Form("%s.EVT%d.DE%d.CLU%d.png",basename,fEventNumber,
-                fDetElemId,fClusterNumber));
+//  if (!fPlot) return;
+//  
+//  TCanvas* c = new TCanvas("MLEM","MLEM",800,600);
+//  c->Draw();
+//  Draw();
+//  c->Modified();
+//  c->Update();
+//  c->Print(Form("%s.EVT%d.DE%d.CLU%d.png",basename,fEventNumber,
+//                fDetElemId,fClusterNumber));
 }
 
 //_____________________________________________________________________________
@@ -764,9 +776,11 @@ AliMUONClusterFinderMLEM::ComputeCoefficients(AliMUONCluster& cluster,
 {
   /// Compute coefficients needed for MLEM algorithm
   
+  Int_t npadTot = cluster.Multiplicity();
   Int_t nPix = fPixArray->GetLast()+1;
   
   //memset(probi,0,nPix*sizeof(Double_t));
+  for (Int_t j = 0; j < npadTot*nPix; ++j) coef[j] = 0.;
   for (Int_t j = 0; j < nPix; ++j) probi[j] = 0.;
 
   Int_t mult = cluster.Multiplicity();
@@ -806,7 +820,7 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
 {
   /// Repeat MLEM algorithm until pixel size becomes sufficiently small
   
-  //  AliCodeTimerAuto("")
+  //  AliCodeTimerAuto("",0)
   
   Int_t nPix = fPixArray->GetLast()+1;
 
@@ -828,7 +842,6 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
     if (cluster.Pad(i)->Status() == fgkZero) ++npadOK;
   }
 
-  TH2D* mlem(0x0);
   Double_t* coef(0x0);
   Double_t* probi(0x0);
   Int_t lc(0); // loop counter
@@ -840,12 +853,12 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
   while (1) 
   {
     ++lc;
-    mlem = (TH2D*) gROOT->FindObject("mlem");
-    delete mlem;
+    delete fHistMlem;
+    fHistMlem = 0x0;
     
     AliDebug(2,Form("lc %d nPix %d(%d) npadTot %d npadOK %d",lc,nPix,fPixArray->GetLast()+1,npadTot,npadOK));
     AliDebug(2,Form("EVT%d PixArray=",fEventNumber));
-    //StdoutToAliDebug(2,fPixArray->Print("","full"));
+    //StdoutToAliDebug(2,fPixArray->Print("full"));
         
     coef = new Double_t [npadTot*nPix];
     probi = new Double_t [nPix];
@@ -895,12 +908,16 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
                     xylim[0],-xylim[1],xylim[2],-xylim[3]
                     ));
     
-    mlem = new TH2D("mlem","mlem",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
+    AliDebug(2,Form("LowestPadCharge=%e",fLowestPadCharge));
+    
+    delete fHistMlem;
+    
+    fHistMlem = new TH2D("mlem","mlem",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
 
     for (Int_t ipix = 0; ipix < nPix; ++ipix) 
     {
-      AliMUONPad* pixPtr = Pixel(ipix);
-      mlem->Fill(pixPtr->Coord(0),pixPtr->Coord(1),pixPtr->Charge());
+      AliMUONPad* pixPtr2 = Pixel(ipix);
+      fHistMlem->Fill(pixPtr2->Coord(0),pixPtr2->Coord(1),pixPtr2->Charge());
     }
 
     // Check if the total charge of pixels is too low
@@ -910,7 +927,7 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
       qTot += Pixel(i)->Charge();
     }
     
-    if ( qTot < 1.e-4 || ( npadOK < 3 && qTot < 7 ) ) 
+    if ( qTot < 1.e-4 || ( npadOK < 3 && qTot < fLowestClusterCharge ) )
     {
       AliDebug(1,Form("Deleting the above cluster (charge %e too low, npadOK=%d)",qTot,npadOK));
       delete [] coef; 
@@ -937,7 +954,7 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
 
     // Calculate position of the center-of-gravity around the maximum pixel
     Double_t xyCOG[2];
-    FindCOG(mlem, xyCOG);
+    FindCOG(xyCOG);
 
     if (TMath::Min(pixPtr->Size(0),pixPtr->Size(1)) < 0.07 && 
         pixPtr->Size(0) > pixPtr->Size(1)) break;
@@ -947,7 +964,7 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
     fPixArray->Sort();
     MaskPeaks(0); // unmask local maxima
     Double_t pixMin = 0.01*Pixel(0)->Charge();
-    pixMin = TMath::Min(pixMin,50.);
+    pixMin = TMath::Min(pixMin,100*fLowestPixelCharge);
 
     // Decrease pixel size and shift pixels to make them centered at 
     // the maximum one
@@ -960,10 +977,10 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
     nPix = 0;
     for (Int_t ipix = 0; ipix < nPix1; ++ipix) 
     {
-      AliMUONPad* pixPtr = Pixel(ipix);
+      AliMUONPad* pixPtr2 = Pixel(ipix);
       if ( nPix >= npadOK  // too many pixels already
            ||
-           pixPtr->Charge() < pixMin && pixPtr->Status() != fgkMustKeep // too low charge
+           ((pixPtr2->Charge() < pixMin) && (pixPtr2->Status() != fgkMustKeep)) // too low charge
            ) 
       { 
         RemovePixel(ipix);
@@ -973,36 +990,36 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
       {
         if (!i) 
         {
-          pixPtr->SetCharge(10);
-          pixPtr->SetSize(indx, pixPtr->Size(indx)/2);
-          width = -pixPtr->Size(indx);
-          pixPtr->Shift(indx, width);
+          pixPtr2->SetCharge(fLowestPadCharge);
+          pixPtr2->SetSize(indx, pixPtr2->Size(indx)/2);
+          width = -pixPtr2->Size(indx);
+          pixPtr2->Shift(indx, width);
           // Shift pixel position
           if (ix) 
           {
             ix = 0;
             for (Int_t j = 0; j < 2; ++j) 
             {
-              shift[j] = pixPtr->Coord(j) - xyCOG[j];
-              shift[j] -= ((Int_t)(shift[j]/pixPtr->Size(j)/2))*pixPtr->Size(j)*2;
+              shift[j] = pixPtr2->Coord(j) - xyCOG[j];
+              shift[j] -= ((Int_t)(shift[j]/pixPtr2->Size(j)/2))*pixPtr2->Size(j)*2;
             }
           } // if (ix)
-          pixPtr->Shift(0, -shift[0]);
-          pixPtr->Shift(1, -shift[1]);
+          pixPtr2->Shift(0, -shift[0]);
+          pixPtr2->Shift(1, -shift[1]);
          ++nPix;
         } 
         else if (nPix < npadOK)
         {
-          pixPtr = new AliMUONPad(*pixPtr);
-          pixPtr->Shift(indx, -2*width);
-         pixPtr->SetStatus(fgkZero);
-          fPixArray->Add(pixPtr);
+          pixPtr2 = new AliMUONPad(*pixPtr2);
+          pixPtr2->Shift(indx, -2*width);
+         pixPtr2->SetStatus(fgkZero);
+          fPixArray->Add(pixPtr2);
          ++nPix;
         } 
        else continue; // skip adjustment of histo limits
         for (Int_t j = 0; j < 4; ++j) 
         {
-          xylim[j] = TMath::Min (xylim[j], (j%2 ? -1 : 1)*pixPtr->Coord(j/2));
+          xylim[j] = TMath::Min (xylim[j], (j%2 ? -1 : 1)*pixPtr2->Coord(j/2));
         }
       } // for (Int_t i=0; i<2;
     } // for (Int_t ipix=0;
@@ -1020,18 +1037,18 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
 
     if (nPix < npadOK)
     {
-      AliMUONPad* pixPtr = Pixel(0);
+      AliMUONPad* pixPtr2 = Pixel(0);
       // add pixels if the maximum is at the limit of pixel area:
       // start from Y-direction
       Int_t j = 0;
       for (Int_t i = 3; i > -1; --i) 
       {
         if (nPix < npadOK && 
-            TMath::Abs((i%2 ? -1 : 1)*xylim[i]-xyCOG[i/2]) < pixPtr->Size(i/2)) 
+            TMath::Abs((i%2 ? -1 : 1)*xylim[i]-xyCOG[i/2]) < pixPtr2->Size(i/2)) 
         {
           //AliMUONPad* p = static_cast<AliMUONPad*>(pixPtr->Clone());
-          AliMUONPad* p = new AliMUONPad(*pixPtr);
-          p->SetCoord(i/2, xyCOG[i/2]+(i%2 ? 2:-2)*pixPtr->Size(i/2));
+          AliMUONPad* p = new AliMUONPad(*pixPtr2);
+          p->SetCoord(i/2, xyCOG[i/2]+(i%2 ? 2:-2)*pixPtr2->Size(i/2));
          xylim[i] = p->Coord(i/2) * (i%2 ? -1 : 1); // update histo limits
           j = TMath::Even (i/2);
           p->SetCoord(j, xyCOG[j]);
@@ -1053,18 +1070,18 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
 
   // remove pixels with low signal or low visibility
   // Cuts are empirical !!!
-  Double_t thresh = TMath::Max (mlem->GetMaximum()/100.,1.);
-  thresh = TMath::Min (thresh,50.);
+  Double_t thresh = TMath::Max (fHistMlem->GetMaximum()/100.,2.0*fLowestPixelCharge);
+  thresh = TMath::Min (thresh,100.0*fLowestPixelCharge);
   Double_t charge = 0;
 
   // Mark pixels which should be removed
   for (Int_t i = 0; i < nPix; ++i) 
   {
-    AliMUONPad* pixPtr = Pixel(i);
-    charge = pixPtr->Charge();
+    AliMUONPad* pixPtr2 = Pixel(i);
+    charge = pixPtr2->Charge();
     if (charge < thresh) 
     {
-      pixPtr->SetCharge(-charge);
+      pixPtr2->SetCharge(-charge);
     }
   }
 
@@ -1072,11 +1089,11 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
   Int_t near = 0;
   for (Int_t i = 0; i < nPix; ++i) 
   {
-    AliMUONPad* pixPtr = Pixel(i);
-    charge = pixPtr->Charge();
+    AliMUONPad* pixPtr2 = Pixel(i);
+    charge = pixPtr2->Charge();
     if (charge > 0) continue;
-    near = FindNearest(pixPtr);
-    pixPtr->SetCharge(0);
+    near = FindNearest(pixPtr2);
+    pixPtr2->SetCharge(0);
     probi[i] = 0; // make it "invisible"
     AliMUONPad* pnear = Pixel(near);
     pnear->SetCharge(pnear->Charge() + (-charge));
@@ -1090,21 +1107,21 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple
   // Update histogram
   for (Int_t i = 0; i < nPix; ++i) 
   {
-    AliMUONPad* pixPtr = Pixel(i);
-    Int_t ix = mlem->GetXaxis()->FindBin(pixPtr->Coord(0));
-    Int_t iy = mlem->GetYaxis()->FindBin(pixPtr->Coord(1));
-    mlem->SetBinContent(ix, iy, pixPtr->Charge());
+    AliMUONPad* pixPtr2 = Pixel(i);
+    Int_t ix = fHistMlem->GetXaxis()->FindBin(pixPtr2->Coord(0));
+    Int_t iy = fHistMlem->GetYaxis()->FindBin(pixPtr2->Coord(1));
+    fHistMlem->SetBinContent(ix, iy, pixPtr2->Charge());
   }
 
   // Try to split into clusters
   Bool_t ok = kTRUE;
-  if (mlem->GetSum() < 1) 
+  if (fHistMlem->GetSum() < 2.0*fLowestPixelCharge) 
   {
     ok = kFALSE;
   }
   else 
   {
-    fSplitter->Split(cluster,mlem,coef,fClusterList);
+    fSplitter->Split(cluster,fHistMlem,coef,fClusterList);
   }
   
   delete [] coef; 
@@ -1131,7 +1148,7 @@ void AliMUONClusterFinderMLEM::MaskPeaks(Int_t mask)
 
 //_____________________________________________________________________________
 void AliMUONClusterFinderMLEM::Mlem(AliMUONCluster& cluster, 
-                                    Double_t* coef, Double_t* probi, 
+                                    const Double_t* coef, Double_t* probi, 
                                     Int_t nIter)
 {
   /// Use MLEM to find pixel charges
@@ -1150,6 +1167,7 @@ void AliMUONClusterFinderMLEM::Mlem(AliMUONCluster& cluster,
     {
       Pixel(ipix)->SetChargeBackup(0);
       // Correct each pixel
+      probi1[ipix] = 0;
       if (probi[ipix] < 0.01) continue; // skip "invisible" pixel
       Double_t sum = 0;
       probi1[ipix] = probMax;
@@ -1176,7 +1194,6 @@ void AliMUONClusterFinderMLEM::Mlem(AliMUONCluster& cluster,
         } 
 
         if (sum1 > 1.e-6) sum += pad->Charge()*coef[indx]/sum1;
-        //if (coef[indx] > 1.e-6) sum += pad->Charge()*coef[indx]/sum1;
       } // for (Int_t j=0;
       AliMUONPad* pixPtr = Pixel(ipix);
       if (probi1[ipix] > 1.e-6) 
@@ -1202,27 +1219,27 @@ void AliMUONClusterFinderMLEM::Mlem(AliMUONCluster& cluster,
 }
 
 //_____________________________________________________________________________
-void AliMUONClusterFinderMLEM::FindCOG(TH2D *mlem, Double_t *xyc)
+void AliMUONClusterFinderMLEM::FindCOG(Double_t *xyc)
 {
   /// Calculate position of the center-of-gravity around the maximum pixel
 
   Int_t ixmax, iymax, ix, nsumx=0, nsumy=0, nsum=0;
   Int_t i1 = -9, j1 = -9;
-  mlem->GetMaximumBin(ixmax,iymax,ix);
-  Int_t nx = mlem->GetNbinsX();
-  Int_t ny = mlem->GetNbinsY();
-  Double_t thresh = mlem->GetMaximum()/10;
+  fHistMlem->GetMaximumBin(ixmax,iymax,ix);
+  Int_t nx = fHistMlem->GetNbinsX();
+  Int_t ny = fHistMlem->GetNbinsY();
+  Double_t thresh = fHistMlem->GetMaximum()/10;
   Double_t x, y, cont, xq=0, yq=0, qq=0;
   
   Int_t ie = TMath::Min(ny,iymax+1), je = TMath::Min(nx,ixmax+1);
   for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
-    y = mlem->GetYaxis()->GetBinCenter(i);
+    y = fHistMlem->GetYaxis()->GetBinCenter(i);
     for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
-      cont = mlem->GetCellContent(j,i);
+      cont = fHistMlem->GetCellContent(j,i);
       if (cont < thresh) continue;
       if (i != i1) {i1 = i; nsumy++;}
       if (j != j1) {j1 = j; nsumx++;}
-      x = mlem->GetXaxis()->GetBinCenter(j);
+      x = fHistMlem->GetXaxis()->GetBinCenter(j);
       xq += x*cont;
       yq += y*cont;
       qq += cont;
@@ -1238,11 +1255,11 @@ void AliMUONClusterFinderMLEM::FindCOG(TH2D *mlem, Double_t *xyc)
     for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
       if (i == iymax) continue;
       for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
-        cont = mlem->GetCellContent(j,i);
+        cont = fHistMlem->GetCellContent(j,i);
         if (cont > cmax) {
           cmax = cont;
-          x = mlem->GetXaxis()->GetBinCenter(j);
-          y = mlem->GetYaxis()->GetBinCenter(i);
+          x = fHistMlem->GetXaxis()->GetBinCenter(j);
+          y = fHistMlem->GetYaxis()->GetBinCenter(i);
           i2 = i;
           j2 = j;
         }
@@ -1262,11 +1279,11 @@ void AliMUONClusterFinderMLEM::FindCOG(TH2D *mlem, Double_t *xyc)
     for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
       if (j == ixmax) continue;
       for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
-        cont = mlem->GetCellContent(j,i);
+        cont = fHistMlem->GetCellContent(j,i);
         if (cont > cmax) {
           cmax = cont;
-          x = mlem->GetXaxis()->GetBinCenter(j);
-          y = mlem->GetYaxis()->GetBinCenter(i);
+          x = fHistMlem->GetXaxis()->GetBinCenter(j);
+          y = fHistMlem->GetYaxis()->GetBinCenter(i);
           i2 = i;
           j2 = j;
         }
@@ -1285,7 +1302,7 @@ void AliMUONClusterFinderMLEM::FindCOG(TH2D *mlem, Double_t *xyc)
 }
 
 //_____________________________________________________________________________
-Int_t AliMUONClusterFinderMLEM::FindNearest(AliMUONPad *pixPtr0)
+Int_t AliMUONClusterFinderMLEM::FindNearest(const AliMUONPad *pixPtr0)
 {
 /// Find the pixel nearest to the given one
 /// (algorithm may be not very efficient)
@@ -1297,7 +1314,7 @@ Int_t AliMUONClusterFinderMLEM::FindNearest(AliMUONPad *pixPtr0)
 
   for (Int_t i = 0; i < nPix; ++i) {
     pixPtr = (AliMUONPad*) fPixArray->UncheckedAt(i);
-    if (pixPtr == pixPtr0 || pixPtr->Charge() < 0.5) continue;
+    if (pixPtr == pixPtr0 || pixPtr->Charge() < fLowestPixelCharge) continue;
     dx = (xc - pixPtr->Coord(0)) / pixPtr->Size(0);
     dy = (yc - pixPtr->Coord(1)) / pixPtr->Size(1);
     r = dx *dx + dy * dy;
@@ -1374,14 +1391,12 @@ Int_t AliMUONClusterFinderMLEM::FindLocalMaxima(TObjArray *pixArray, Int_t *loca
 
   AliDebug(1,Form("nPix=%d",pixArray->GetLast()+1));
 
-  TH2D *hist = NULL;
-  //if (pixArray == fPixArray) hist = (TH2D*) gROOT->FindObject("anode");
-  //else { hist = (TH2D*) gROOT->FindObject("anode1"); cout << hist << endl; }
-  //if (hist) hist->Delete();
   Double_t xylim[4] = {999, 999, 999, 999};
 
   Int_t nPix = pixArray->GetEntriesFast();
+  
+  if ( nPix <= 0 ) return 0;
+  
   AliMUONPad *pixPtr = 0;
   for (Int_t ipix = 0; ipix < nPix; ++ipix) {
     pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix);
@@ -1392,11 +1407,11 @@ Int_t AliMUONClusterFinderMLEM::FindLocalMaxima(TObjArray *pixArray, Int_t *loca
 
   Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
   Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
-  if (pixArray == fPixArray) hist = new TH2D("anode","anode",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
-  else hist = new TH2D("anode1","anode1",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
+  if (pixArray == fPixArray) fHistAnode = new TH2D("anode","anode",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
+  else fHistAnode = new TH2D("anode1","anode1",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
   for (Int_t ipix = 0; ipix < nPix; ++ipix) {
     pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix);
-    hist->Fill(pixPtr->Coord(0), pixPtr->Coord(1), pixPtr->Charge());
+    fHistAnode->Fill(pixPtr->Coord(0), pixPtr->Coord(1), pixPtr->Charge());
   }
 //  if (fDraw && pixArray == fPixArray) fDraw->DrawHist("c2", hist);
 
@@ -1407,10 +1422,10 @@ Int_t AliMUONClusterFinderMLEM::FindLocalMaxima(TObjArray *pixArray, Int_t *loca
   for (Int_t i = 1; i <= ny; ++i) {
     indx = (i-1) * nx;
     for (Int_t j = 1; j <= nx; ++j) {
-      if (hist->GetCellContent(j,i) < 0.5) continue;
+      if (fHistAnode->GetCellContent(j,i) < fLowestPixelCharge) continue;
       //if (isLocalMax[indx+j-1] < 0) continue;
       if (isLocalMax[indx+j-1] != 0) continue;
-      FlagLocalMax(hist, i, j, isLocalMax);
+      FlagLocalMax(fHistAnode, i, j, isLocalMax);
     }
   }
 
@@ -1419,11 +1434,15 @@ Int_t AliMUONClusterFinderMLEM::FindLocalMaxima(TObjArray *pixArray, Int_t *loca
     for (Int_t j = 1; j <= nx; ++j) {
       if (isLocalMax[indx+j-1] > 0) { 
        localMax[nMax] = indx + j - 1; 
-       maxVal[nMax++] = hist->GetCellContent(j,i);
-       ((AliMUONPad*)fSplitter->BinToPix(hist, j, i))->SetStatus(fgkMustKeep);
-       if (nMax > 99) AliFatal(" Too many local maxima !!!");
+       maxVal[nMax++] = fHistAnode->GetCellContent(j,i);
+       ((AliMUONPad*)fSplitter->BinToPix(fHistAnode, j, i))->SetStatus(fgkMustKeep);
+       if (nMax > 99) break;
       }
     }
+    if (nMax > 99) {
+      AliError(" Too many local maxima !!!");
+      break;
+    }
   }
   if (fDebug) cout << " Local max: " << nMax << endl;
   delete [] isLocalMax; 
@@ -1466,12 +1485,11 @@ void AliMUONClusterFinderMLEM::FlagLocalMax(TH2D *hist, Int_t i, Int_t j, Int_t
 
 //_____________________________________________________________________________
 void AliMUONClusterFinderMLEM::FindCluster(AliMUONCluster& cluster, 
-                                           Int_t *localMax, Int_t iMax)
+                                           const Int_t *localMax, Int_t iMax)
 {
 /// Find pixel cluster around local maximum \a iMax and pick up pads
 /// overlapping with it
 
-  TH2D *hist = (TH2D*) gROOT->FindObject("anode");
   /* Just for check
   TCanvas* c = new TCanvas("Anode","Anode",800,600);
   c->cd();
@@ -1480,8 +1498,8 @@ void AliMUONClusterFinderMLEM::FindCluster(AliMUONCluster& cluster,
   Int_t tmp;
   cin >> tmp;
   */
-  Int_t nx = hist->GetNbinsX();
-  Int_t ny = hist->GetNbinsY();
+  Int_t nx = fHistAnode->GetNbinsX();
+  Int_t ny = fHistAnode->GetNbinsY();
   Int_t ic = localMax[iMax] / nx + 1;
   Int_t jc = localMax[iMax] % nx + 1;
   Int_t nxy = ny * nx;
@@ -1491,14 +1509,14 @@ void AliMUONClusterFinderMLEM::FindCluster(AliMUONCluster& cluster,
   // Drop all pixels from the array - pick up only the ones from the cluster
   fPixArray->Delete();
 
-  Double_t wx = hist->GetXaxis()->GetBinWidth(1)/2; 
-  Double_t wy = hist->GetYaxis()->GetBinWidth(1)/2;  
-  Double_t yc = hist->GetYaxis()->GetBinCenter(ic);
-  Double_t xc = hist->GetXaxis()->GetBinCenter(jc);
-  Double_t cont = hist->GetCellContent(jc,ic);
+  Double_t wx = fHistAnode->GetXaxis()->GetBinWidth(1)/2; 
+  Double_t wy = fHistAnode->GetYaxis()->GetBinWidth(1)/2;  
+  Double_t yc = fHistAnode->GetYaxis()->GetBinCenter(ic);
+  Double_t xc = fHistAnode->GetXaxis()->GetBinCenter(jc);
+  Double_t cont = fHistAnode->GetCellContent(jc,ic);
   fPixArray->Add(new AliMUONPad (xc, yc, wx, wy, cont));
   used[(ic-1)*nx+jc-1] = kTRUE;
-  AddBinSimple(hist, ic, jc);
+  AddBinSimple(fHistAnode, ic, jc);
   //fSplitter->AddBin(hist, ic, jc, 1, used, (TObjArray*)0); // recursive call
 
   Int_t nPix = fPixArray->GetEntriesFast();
@@ -1540,23 +1558,23 @@ void AliMUONClusterFinderMLEM::FindCluster(AliMUONCluster& cluster,
 
 //_____________________________________________________________________________
 void 
-AliMUONClusterFinderMLEM::AddBinSimple(TH2D *mlem, Int_t ic, Int_t jc)
+AliMUONClusterFinderMLEM::AddBinSimple(TH2D *hist, Int_t ic, Int_t jc)
 {
   /// Add adjacent bins (+-1 in X and Y) to the cluster
   
-  Int_t nx = mlem->GetNbinsX();
-  Int_t ny = mlem->GetNbinsY();
-  Double_t cont1, cont = mlem->GetCellContent(jc,ic);
+  Int_t nx = hist->GetNbinsX();
+  Int_t ny = hist->GetNbinsY();
+  Double_t cont1, cont = hist->GetCellContent(jc,ic);
   AliMUONPad *pixPtr = 0;
   
   Int_t ie = TMath::Min(ic+1,ny), je = TMath::Min(jc+1,nx);
   for (Int_t i = TMath::Max(ic-1,1); i <= ie; ++i) {
     for (Int_t j = TMath::Max(jc-1,1); j <= je; ++j) {
-      cont1 = mlem->GetCellContent(j,i);
+      cont1 = hist->GetCellContent(j,i);
       if (cont1 > cont) continue;
-      if (cont1 < 0.5) continue;
-      pixPtr = new AliMUONPad (mlem->GetXaxis()->GetBinCenter(j), 
-                              mlem->GetYaxis()->GetBinCenter(i), 0, 0, cont1);
+      if (cont1 < fLowestPixelCharge) continue;
+      pixPtr = new AliMUONPad (hist->GetXaxis()->GetBinCenter(j), 
+                              hist->GetYaxis()->GetBinCenter(i), 0, 0, cont1);
       fPixArray->Add(pixPtr);
     }
   }
@@ -1594,19 +1612,19 @@ void AliMUONClusterFinderMLEM::AddVirtualPad(AliMUONCluster& cluster)
   Bool_t same = kFALSE;
   if (TMath::Abs(dim0.Y()-dim1.Y()) < fgkDistancePrecision) same = kTRUE; // the same pad size on both planes 
 
-  AliMpIntPair cn;
+  Long_t cn;
   Bool_t check[2] = {kFALSE, kFALSE};
   Int_t nxy[2];
   nxy[0] = nxy[1] = 0;
   for (Int_t inb = 0; inb < 2; ++inb) {
     cn = cluster.NofPads(nonb[inb], 0, kTRUE);
-    if (inb == 0 && cn.GetFirst() == 2) check[inb] = kTRUE; // check non-bending plane
-    else if (inb == 1 && cn.GetSecond() == 2) check[inb] = kTRUE; // check bending plane
+    if (inb == 0 && AliMp::PairFirst(cn) == 2) check[inb] = kTRUE; // check non-bending plane
+    else if (inb == 1 && AliMp::PairSecond(cn) == 2) check[inb] = kTRUE; // check bending plane
     if (same) {
-      nxy[0] = TMath::Max (nxy[0], cn.GetFirst());
-      nxy[1] = TMath::Max (nxy[1], cn.GetSecond());
+      nxy[0] = TMath::Max (nxy[0], AliMp::PairFirst(cn));
+      nxy[1] = TMath::Max (nxy[1], AliMp::PairSecond(cn));
       if (inb == 0 && nxy[0] < 2) nonb[inb] = !nonb[inb];
-      else if (inb == 1 && cn.GetSecond() < 2) nonb[inb] = !nonb[inb];
+      else if (inb == 1 && AliMp::PairSecond(cn) < 2) nonb[inb] = !nonb[inb];
     }
   }
   if (same) {
@@ -1668,17 +1686,16 @@ void AliMUONClusterFinderMLEM::AddVirtualPad(AliMUONCluster& cluster)
       TVector2 pos;
       if (inb == 0) pos.Set (pad->X() + idir * (pad->DX()+fgkDistancePrecision), pad->Y());
       else pos.Set (pad->X(), pad->Y() + idir * (pad->DY()+fgkDistancePrecision));
-      //AliMpPad mppad = fSegmentation[nonb[inb]]->PadByPosition(pos,kTRUE);
-      AliMpPad mppad = fSegmentation[nonb[inb]]->PadByPosition(pos,kFALSE);
+      //AliMpPad mppad = fkSegmentation[nonb[inb]]->PadByPosition(pos,kTRUE);
+      AliMpPad mppad = fkSegmentation[nonb[inb]]->PadByPosition(pos.X(), pos.Y(),kFALSE);
       if (!mppad.IsValid()) continue; // non-existing pad
-      cn = mppad.GetIndices();
-      AliMUONPad muonPad(fDetElemId, nonb[inb], cn.GetFirst(), cn.GetSecond(), 
-                        mppad.Position().X(), mppad.Position().Y(), 
-                        mppad.Dimensions().X(), mppad.Dimensions().Y(), 0);
-      if (inb == 0) muonPad.SetCharge(TMath::Min (amax[j]/100, 5.));
+      AliMUONPad muonPad(fDetElemId, nonb[inb], mppad.GetIx(), mppad.GetIy(), 
+                        mppad.GetPositionX(), mppad.GetPositionY(), 
+                        mppad.GetDimensionX(), mppad.GetDimensionY(), 0);
+      if (inb == 0) muonPad.SetCharge(TMath::Min (amax[j]/100, fLowestPadCharge));
       //else muonPad.SetCharge(TMath::Min (amax[j]/15, fgkZeroSuppression));
-      else muonPad.SetCharge(TMath::Min (amax[j]/15, 6.));
-      if (muonPad.Charge() < 1.) muonPad.SetCharge(1.);
+      else muonPad.SetCharge(TMath::Min (amax[j]/15, fLowestPadCharge));
+      if (muonPad.Charge() < 2.0*fLowestPixelCharge) muonPad.SetCharge(2.0*fLowestPixelCharge);
       muonPad.SetReal(kFALSE);
       if (fDebug) printf(" ***** Add virtual pad in %d direction ***** %f %f %f %3d %3d %f %f \n",
                         inb, muonPad.Charge(), muonPad.X(), muonPad.Y(), muonPad.Ix(), 
@@ -1705,10 +1722,10 @@ void AliMUONClusterFinderMLEM::PadsInXandY(AliMUONCluster& cluster,
        
   Bool_t mustMatch(kTRUE);
 
-  AliMpIntPair cn = cluster.NofPads(statusToTest,mustMatch);
+  Long_t cn = cluster.NofPads(statusToTest,mustMatch);
   
-  nInX = cn.GetFirst();
-  nInY = cn.GetSecond();
+  nInX = AliMp::PairFirst(cn);
+  nInY = AliMp::PairSecond(cn);
 }
 
 //_____________________________________________________________________________
@@ -1748,7 +1765,7 @@ void AliMUONClusterFinderMLEM::Simple(AliMUONCluster& cluster)
     */
     if (!pad->IsSaturated()) pad->SetStatus(fgkUseForFit);
   }
-  nfit = fSplitter->Fit(cluster,1, nForFit, clustFit, clusters, parOk, fClusterList);
+  nfit = fSplitter->Fit(cluster,1, nForFit, clustFit, clusters, parOk, fClusterList, fHistMlem);
 }
 
 //_____________________________________________________________________________
@@ -1772,4 +1789,14 @@ AliMUONClusterFinderMLEM::Print(Option_t* what) const
   }
 }
 
+//_____________________________________________________________________________
+void 
+AliMUONClusterFinderMLEM::SetChargeHints(Double_t lowestPadCharge, Double_t lowestClusterCharge)
+{
+  /// Set some thresholds we use later on in the algorithm
+  fLowestPadCharge=lowestPadCharge;
+  fLowestClusterCharge=lowestClusterCharge;
+  fLowestPixelCharge=fLowestPadCharge/12.0; 
+}
+