]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONClusterFinderVS.cxx
Stored additional energy information
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderVS.cxx
index 2b95fb300fc14a8474fd88d21e063bf4f5cb90d1..0c5f92825628fb0f69e09da260413b8dce934e50 100644 (file)
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-/*
-$Log$
-Revision 1.11  2000/10/06 09:04:05  morsch
-- Dummy z-arguments in GetPadI, SetHit, FirstPad replaced by real z-coordinate
-       to make code work with slat chambers (AM)
-- Replace GetPadI calls with unchecked x,y coordinates by pad iterator calls wherever possible.
-
-Revision 1.10  2000/10/03 13:51:57  egangler
-Removal of useless dependencies via forward declarations
-
-Revision 1.9  2000/10/02 16:58:29  egangler
-Cleaning of the code :
--> coding conventions
--> void Streamers
--> some useless includes removed or replaced by "class" statement
-
-Revision 1.8  2000/07/03 11:54:57  morsch
-AliMUONSegmentation and AliMUONHitMap have been replaced by AliSegmentation and AliHitMap in STEER
-The methods GetPadIxy and GetPadXxy of AliMUONSegmentation have changed name to GetPadI and GetPadC.
-
-Revision 1.7  2000/06/28 15:16:35  morsch
-(1) Client code adapted to new method signatures in AliMUONSegmentation (see comments there)
-to allow development of slat-muon chamber simulation and reconstruction code in the MUON
-framework. The changes should have no side effects (mostly dummy arguments).
-(2) Hit disintegration uses 3-dim hit coordinates to allow simulation
-of chambers with overlapping modules (MakePadHits, Disintegration).
-
-Revision 1.6  2000/06/28 12:19:18  morsch
-More consequent seperation of global input data services (AliMUONClusterInput singleton) and the
-cluster and hit reconstruction algorithms in AliMUONClusterFinderVS.
-AliMUONClusterFinderVS becomes the base class for clustering and hit reconstruction.
-It requires two cathode planes. Small modifications in the code will make it usable for
-one cathode plane and, hence, more general (for test beam data).
-AliMUONClusterFinder is now obsolete.
-
-Revision 1.5  2000/06/28 08:06:10  morsch
-Avoid global variables in AliMUONClusterFinderVS by seperating the input data for the fit from the
-algorithmic part of the class. Input data resides inside the AliMUONClusterInput singleton.
-It also naturally takes care of the TMinuit instance.
-
-Revision 1.4  2000/06/27 16:18:47  gosset
-Finally correct implementation of xm, ym, ixm, iym sizes
-when at least three local maxima on cathode 1 or on cathode 2
-
-Revision 1.3  2000/06/22 14:02:45  morsch
-Parameterised size of xm[], ym[], ixm[], iym[] correctly implemented (PH)
-Some HP scope problems corrected (PH)
-
-Revision 1.2  2000/06/15 07:58:48  morsch
-Code from MUON-dev joined
-
-Revision 1.1.2.3  2000/06/09 21:58:33  morsch
-Most coding rule violations corrected.
-
-Revision 1.1.2.2  2000/02/15 08:33:52  morsch
-Error in calculation of contribution map for double clusters (Split method) corrected   (A.M.)
-Error in determination of track list for double cluster (FillCluster method) corrected  (A.M.)
-Revised and extended SplitByLocalMaxima method (Isabelle Chevrot):
-       - For clusters with more than 2 maxima on one of the cathode planes all valid
-       combinations of maxima on the two cathodes are preserved. The position of the maxima is
-       taken as the hit position.
-       - New FillCluster method with 2 arguments to find tracks associated to the clusters
-       defined above added. (Method destinction by argument list not very elegant in this case,
-       should be revides (A.M.)
-       - Bug in if-statement to handle maximum 1 maximum per plane corrected
-       - Two cluster per cathode but only 1 combination valid is handled.
-       - More rigerous treatment of 1-2 and 2-1 combinations of maxima.
 
-*/
+/* $Id$ */
 
 #include "AliMUONClusterFinderVS.h"
 #include "AliMUONDigit.h"
@@ -102,7 +35,7 @@ Revised and extended SplitByLocalMaxima method (Isabelle Chevrot):
 #include <TF1.h>
 
 #include <stdio.h>
-#include <iostream.h>
+#include <Riostream.h>
 
 //_____________________________________________________________________
 // This function is minimized in the double-Mathieson fit
@@ -120,10 +53,18 @@ ClassImp(AliMUONClusterFinderVS)
     fHitMap[0] = 0;
     fHitMap[1] = 0;
     fTrack[0]=fTrack[1]=-1;
+    fDebugLevel = 0; // make silent default
+    fGhostChi2Cut = 1e6; // nothing done by default
+    fSeg[0]    = 0;
+    fSeg[1]    = 0;
+    for(Int_t i=0; i<100; i++) {
+      for (Int_t j=0; j<2; j++) {
+        fDig[i][j] = 0;
+      }
+    }
 }
 
-AliMUONClusterFinderVS::AliMUONClusterFinderVS(
-    const AliMUONClusterFinderVS & clusterFinder)
+AliMUONClusterFinderVS::AliMUONClusterFinderVS(const AliMUONClusterFinderVS & clusterFinder):TObject(clusterFinder)
 {
 // Dummy copy Constructor
     ;
@@ -158,10 +99,10 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
            // pointer to digit
            fDig[i][cath]=fInput->Digit(cath, c->fIndexMap[i][cath]);
            // pad coordinates
-           fIx[i][cath]= fDig[i][cath]->fPadX;
-           fIy[i][cath]= fDig[i][cath]->fPadY;
+           fIx[i][cath]= fDig[i][cath]->PadX();
+           fIy[i][cath]= fDig[i][cath]->PadY();
            // pad charge
-           fQ[i][cath] = fDig[i][cath]->fSignal;
+           fQ[i][cath] = fDig[i][cath]->Signal();
            // pad centre coordinates
            fSeg[cath]->
                GetPadC(fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath], fZ[i][cath]);
@@ -179,7 +120,6 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 //  +++++++++++++++++++++++++++++++*************++++++++
     if ((fNLocal[0]==1 && (fNLocal[1]==0 ||  fNLocal[1]==1)) || 
        (fNLocal[0]==0 && fNLocal[1]==1)) {
-
 // Perform combined single Mathieson fit
 // Initial values for coordinates (x,y) 
 
@@ -196,14 +136,16 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
            fXInit[0]=c->fX[1];
            fYInit[0]=c->fY[1];
        }
-       fprintf(stderr,"\n cas (1) CombiSingleMathiesonFit(c)\n");
+       if (fDebugLevel)
+           fprintf(stderr,"\n cas (1) CombiSingleMathiesonFit(c)\n");
        chi2=CombiSingleMathiesonFit(c);
 //     Int_t ndf = fgNbins[0]+fgNbins[1]-2;
 //     Float_t prob = TMath::Prob(Double_t(chi2),ndf);
 //     prob1->Fill(prob);
 //     chi2_1->Fill(chi2);
        oldchi2=chi2;
-       fprintf(stderr," chi2 %f ",chi2);
+       if (fDebugLevel)
+           fprintf(stderr," chi2 %f ",chi2);
 
        c->fX[0]=fXFit[0];
        c->fY[0]=fYFit[0];
@@ -212,11 +154,11 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
        c->fY[1]=fYFit[0];
        c->fChi2[0]=chi2;
        c->fChi2[1]=chi2;
+        // Force on anod
        c->fX[0]=fSeg[0]->GetAnod(c->fX[0]);
        c->fX[1]=fSeg[1]->GetAnod(c->fX[1]);
        
 // If reasonable chi^2 add result to the list of rawclusters
-       //      if (chi2 < 50) {
        if (chi2 < 0.3) {
            AddRawCluster(*c);
 // If not try combined double Mathieson Fit
@@ -242,6 +184,7 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 //  Initial value for charge ratios
            fQrInit[0]=0.5;
            fQrInit[1]=0.5;
+           if (fDebugLevel)
            fprintf(stderr,"\n cas (1) CombiDoubleMathiesonFit(c)\n");
            chi2=CombiDoubleMathiesonFit(c);
 //         Int_t ndf = fgNbins[0]+fgNbins[1]-6;
@@ -268,7 +211,7 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
     } else if (fNLocal[0]==2 &&  fNLocal[1]==2) {
 //
 //  Let's look for ghosts first 
-//
+
        Float_t xm[4][2], ym[4][2];
        Float_t dpx, dpy, dx, dy;
        Int_t ixm[4][2], iym[4][2];
@@ -298,10 +241,15 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 
 // Analyse the combinations and keep those that are possible !
 // For each combination check consistency in x and y   
-       Int_t iacc;
-       Bool_t accepted[4];
+       Int_t   iacc;
+       Bool_t  accepted[4];
+       Float_t dr[4] = {1.e4, 1.e4, 1.e4, 1.e4};
        iacc=0;
-       
+
+// In case of staggering maxima are displaced by exactly half the pad-size in y. 
+// We have to take into account the numerical precision in the consistency check;      
+       Float_t eps = 1.e-5;
+//
        for (ico=0; ico<4; ico++) {
            accepted[ico]=kFALSE;
 // cathode one: x-coordinate
@@ -312,23 +260,62 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
            isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
            dpy=fSeg[1]->Dpy(isec)/2.;
            dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
-//         printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
-           if ((dx <= dpx) && (dy <= dpy)) {
+           if (fDebugLevel>1) 
+               printf("\n %i %f %f %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy, dx, dpx );
+           if ((dx <= dpx) && (dy <= dpy+eps)) {
                // consistent
                accepted[ico]=kTRUE;
+               dr[ico] = TMath::Sqrt(dx*dx+dy*dy);
                iacc++;
            } else {
                // reject
                accepted[ico]=kFALSE;
            }
        }
+       printf("\n iacc= %d:\n", iacc);
+       if (iacc == 3) {
+           if (accepted[0] && accepted[1]) {
+               if (dr[0] >= dr[1]) {
+                   accepted[0]=kFALSE;
+               } else {
+                   accepted[1]=kFALSE;
+               }
+           }
 
-       if (iacc==2) {
-           fprintf(stderr,"\n iacc=2: No problem ! \n");
-       } else if (iacc==4) {
-           fprintf(stderr,"\n iacc=4: Ok, but ghost problem !!! \n");
-       } else if (iacc==0) {
-           fprintf(stderr,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
+           if (accepted[2] && accepted[3]) {
+               if (dr[2] >= dr[3]) {
+                   accepted[2]=kFALSE;
+               } else {
+                   accepted[3]=kFALSE;
+               }
+           }
+/*         
+// eliminate one candidate
+           Float_t drmax = 0;
+           Int_t icobad = -1;
+
+           for (ico=0; ico<4; ico++) {
+               if (accepted[ico] && dr[ico] > drmax) {
+                   icobad = ico;
+                   drmax  = dr[ico];
+               }
+           }
+           
+           accepted[icobad] = kFALSE;
+*/
+           iacc = 2;
+       }
+       
+       
+       printf("\n iacc= %d:\n", iacc);
+       if (fDebugLevel) {
+           if (iacc==2) {
+               fprintf(stderr,"\n iacc=2: No problem ! \n");
+           } else if (iacc==4) {
+               fprintf(stderr,"\n iacc=4: Ok, but ghost problem !!! \n");
+           } else if (iacc==0) {
+               fprintf(stderr,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
+           }
        }
 
 //  Initial value for charge ratios
@@ -347,7 +334,6 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 // ******* iacc = 1 *******
 // Only one combination found between the 2 cathodes
        if (iacc==1) {
-
 // Initial values for the 2 maxima (x,y)
 
 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
@@ -377,13 +363,15 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
                fXInit[1]=xm[0][0];
                fYInit[1]=ym[0][0];
            }
-           fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
+           if (fDebugLevel)
+               fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
            chi2=CombiDoubleMathiesonFit(c);
 //         Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 //         Float_t prob = TMath::Prob(chi2,ndf);
 //         prob2->Fill(prob);
 //         chi2_2->Fill(chi2);
-           fprintf(stderr," chi2 %f\n",chi2);
+           if (fDebugLevel)
+               fprintf(stderr," chi2 %f\n",chi2);
 
 // If reasonable chi^2 add result to the list of rawclusters
            if (chi2<10) {
@@ -417,13 +405,15 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
                    fXInit[1]=xm[0][1];
                    fYInit[1]=ym[0][1];
                }
-               fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
+               if (fDebugLevel)
+                   fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
                chi2=CombiDoubleMathiesonFit(c);
 //             Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 //             Float_t prob = TMath::Prob(chi2,ndf);
 //             prob2->Fill(prob);
 //             chi2_2->Fill(chi2);
-               fprintf(stderr," chi2 %f\n",chi2);
+               if (fDebugLevel)
+                   fprintf(stderr," chi2 %f\n",chi2);
 
 // If reasonable chi^2 add result to the list of rawclusters
                if (chi2<10) {
@@ -437,6 +427,8 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
                            for (cath=0; cath<2; cath++) {
                                cnew.fX[cath]=Float_t(xm[ico][1]);
                                cnew.fY[cath]=Float_t(ym[ico][0]);
+                               cnew.fZ[cath]=fZPlane;
+                               
                                cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
                                for (i=0; i<fMul[cath]; i++) {
                                    cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
@@ -458,12 +450,11 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 // ******* iacc = 2 *******
 // Two combinations found between the 2 cathodes
        if (iacc==2) {
-
 // Was the same maximum taken twice
            if ((accepted[0]&&accepted[1]) || (accepted[2]&&accepted[3])) {
                fprintf(stderr,"\n Maximum taken twice !!!\n");
 
-// Have a try !! with that 
+// Have a try !! with that
                if (accepted[0]&&accepted[3]) {
                    fXInit[0]=xm[0][1];
                    fYInit[0]=ym[0][0];
@@ -475,7 +466,8 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
                    fXInit[1]=xm[3][1];
                    fYInit[1]=ym[3][0];
                }
-               fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
+               if (fDebugLevel)
+                   fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
                chi2=CombiDoubleMathiesonFit(c);
 //                 Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 //                 Float_t prob = TMath::Prob(chi2,ndf);
@@ -496,13 +488,15 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
                    fXInit[1]=xm[2][1];
                    fYInit[1]=ym[2][0];
                }
-               fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
+               if (fDebugLevel)
+                   fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
                chi2=CombiDoubleMathiesonFit(c);
 //                 Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 //                 Float_t prob = TMath::Prob(chi2,ndf);
 //                 prob2->Fill(prob);
 //                 chi2_2->Fill(chi2);
-               fprintf(stderr," chi2 %f\n",chi2);
+               if (fDebugLevel)
+                   fprintf(stderr," chi2 %f\n",chi2);
                Split(c);
            }
            
@@ -511,30 +505,90 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 // Ghost !!
        } else if (iacc==4) {
 // Perform fits for the two possibilities !!   
+// Accept if charges are compatible on both cathodes
+// If none are compatible, keep everything
            fXInit[0]=xm[0][1];
            fYInit[0]=ym[0][0];
            fXInit[1]=xm[3][1];
            fYInit[1]=ym[3][0];
-           fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
+           if (fDebugLevel)
+               fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
            chi2=CombiDoubleMathiesonFit(c);
 //             Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 //             Float_t prob = TMath::Prob(chi2,ndf);
 //             prob2->Fill(prob);
 //             chi2_2->Fill(chi2);
-           fprintf(stderr," chi2 %f\n",chi2);
-           Split(c);
+           if (fDebugLevel)
+               fprintf(stderr," chi2 %f\n",chi2);
+           // store results of fit and postpone decision
+           Double_t sXFit[2],sYFit[2],sQrFit[2];
+           Float_t sChi2[2];
+           for (Int_t i=0;i<2;i++) {
+               sXFit[i]=fXFit[i];
+               sYFit[i]=fYFit[i];
+               sQrFit[i]=fQrFit[i];
+               sChi2[i]=fChi2[i];
+           }
            fXInit[0]=xm[1][1];
            fYInit[0]=ym[1][0];
            fXInit[1]=xm[2][1];
            fYInit[1]=ym[2][0];
-           fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
+           if (fDebugLevel)
+               fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
            chi2=CombiDoubleMathiesonFit(c);
 //             ndf = fgNbins[0]+fgNbins[1]-6;
 //             prob = TMath::Prob(chi2,ndf);
 //             prob2->Fill(prob);
 //             chi2_2->Fill(chi2);
-           fprintf(stderr," chi2 %f\n",chi2);
-           Split(c);
+           if (fDebugLevel)
+               fprintf(stderr," chi2 %f\n",chi2);
+           // We have all informations to perform the decision
+           // Compute the chi2 for the 2 possibilities
+           Float_t chi2fi,chi2si,chi2f,chi2s;
+
+           chi2f = (TMath::Log(fInput->TotalCharge(0)*fQrFit[0]
+                 /  (fInput->TotalCharge(1)*fQrFit[1]) )
+                 / fInput->Response()->ChargeCorrel() );
+           chi2f *=chi2f;
+           chi2fi = (TMath::Log(fInput->TotalCharge(0)*(1-fQrFit[0])
+                 /  (fInput->TotalCharge(1)*(1-fQrFit[1])) )
+                 / fInput->Response()->ChargeCorrel() );
+           chi2f += chi2fi*chi2fi;
+
+           chi2s = (TMath::Log(fInput->TotalCharge(0)*sQrFit[0]
+                 /  (fInput->TotalCharge(1)*sQrFit[1]) )
+                 / fInput->Response()->ChargeCorrel() );
+           chi2s *=chi2s;
+           chi2si = (TMath::Log(fInput->TotalCharge(0)*(1-sQrFit[0])
+                 /  (fInput->TotalCharge(1)*(1-sQrFit[1])) )
+                 / fInput->Response()->ChargeCorrel() );
+           chi2s += chi2si*chi2si;
+
+           // usefull to store the charge matching chi2 in the cluster
+           // fChi2[0]=sChi2[1]=chi2f;
+           // fChi2[1]=sChi2[0]=chi2s;
+
+           if (chi2f<=fGhostChi2Cut && chi2s<=fGhostChi2Cut)
+               c->fGhost=1;
+           if   (chi2f>fGhostChi2Cut && chi2s>fGhostChi2Cut) {
+               // we keep the ghost
+               c->fGhost=2;
+               chi2s=-1;
+               chi2f=-1;
+           }
+           if (chi2f<=fGhostChi2Cut)
+               Split(c);
+           if (chi2s<=fGhostChi2Cut) {
+               // retreive saved values
+               for (Int_t i=0;i<2;i++) {
+                   fXFit[i]=sXFit[i];
+                   fYFit[i]=sYFit[i];
+                   fQrFit[i]=sQrFit[i];
+                   fChi2[i]=sChi2[i];
+               }
+               Split(c);
+           }
+           c->fGhost=0;
        }
 
     } else if (fNLocal[0]==2 &&  fNLocal[1]==1) {
@@ -570,7 +624,11 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
        Int_t iacc;
        Bool_t accepted[4];
        iacc=0;
+       // In case of staggering maxima are displaced by exactly half the pad-size in y. 
+        // We have to take into account the numerical precision in the consistency check;
        
+       Float_t eps = 1.e-5;
+
        for (ico=0; ico<2; ico++) {
            accepted[ico]=kFALSE;
            isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
@@ -579,8 +637,9 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
            isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
            dpy=fSeg[1]->Dpy(isec)/2.;
            dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
-//         printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
-           if ((dx <= dpx) && (dy <= dpy)) {
+           if (fDebugLevel>1)
+               printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
+           if ((dx <= dpx) && (dy <= dpy+eps)) {
                // consistent
                accepted[ico]=kTRUE;
                iacc++;
@@ -592,8 +651,31 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
        
        Float_t chi21 = 100;
        Float_t chi22 = 100;
+       Float_t chi23 = 100;
+
+       //  Initial value for charge ratios
+       fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
+           Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
+       fQrInit[1]=fQrInit[0];
        
-       if (accepted[0]) {
+       if (accepted[0] && accepted[1]) {
+           
+           fXInit[0]=0.5*(xm[0][1]+xm[0][0]);
+           fYInit[0]=ym[0][0];
+           fXInit[1]=0.5*(xm[0][1]+xm[1][0]);
+           fYInit[1]=ym[1][0];
+           fQrInit[0]=0.5;
+           fQrInit[1]=0.5;
+           chi23=CombiDoubleMathiesonFit(c);
+           if (chi23<10) {
+               Split(c);
+               Float_t yst;
+               yst = fYFit[0];
+               fYFit[0] = fYFit[1];
+               fYFit[1] = yst;
+               Split(c);
+           }
+       } else if (accepted[0]) {
            fXInit[0]=xm[0][1];
            fYInit[0]=ym[0][0];
            fXInit[1]=xm[1][0];
@@ -603,7 +685,8 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 //         Float_t prob = TMath::Prob(chi2,ndf);
 //         prob2->Fill(prob);
 //         chi2_2->Fill(chi21);
-           fprintf(stderr," chi2 %f\n",chi21);
+           if (fDebugLevel)
+               fprintf(stderr," chi2 %f\n",chi21);
            if (chi21<10) Split(c);
        } else if (accepted[1]) {
            fXInit[0]=xm[1][1];
@@ -615,11 +698,12 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 //         Float_t prob = TMath::Prob(chi2,ndf);
 //         prob2->Fill(prob);
 //         chi2_2->Fill(chi22);
-           fprintf(stderr," chi2 %f\n",chi22);
+           if (fDebugLevel)
+               fprintf(stderr," chi2 %f\n",chi22);
            if (chi22<10) Split(c);
        }
 
-       if (chi21 > 10 && chi22 > 10) {
+       if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
 // We keep only the combination found (X->cathode 2, Y->cathode 1)
            for (Int_t ico=0; ico<2; ico++) {
                if (accepted[ico]) {
@@ -628,6 +712,7 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
                    for (cath=0; cath<2; cath++) {
                        cnew.fX[cath]=Float_t(xm[ico][1]);
                        cnew.fY[cath]=Float_t(ym[ico][0]);
+                       cnew.fZ[cath]=fZPlane;
                        cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
                        for (i=0; i<fMul[cath]; i++) {
                            cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
@@ -648,7 +733,6 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 //  (3') One local maximum on cathode 1 and two maxima on cathode 2 
 //  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     } else if (fNLocal[0]==1 && fNLocal[1]==2) {
-       
        Float_t xm[4][2], ym[4][2];
        Float_t dpx, dpy, dx, dy;
        Int_t ixm[4][2], iym[4][2];
@@ -677,6 +761,10 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
        Int_t iacc;
        Bool_t accepted[4];
        iacc=0;
+        // In case of staggering maxima are displaced by exactly half the pad-size in y. 
+        // We have to take into account the numerical precision in the consistency check;      
+       Float_t eps = 1.e-5;
+
        
        for (ico=0; ico<2; ico++) {
            accepted[ico]=kFALSE;
@@ -686,8 +774,9 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
            isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
            dpy=fSeg[1]->Dpy(isec)/2.;
            dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
-//         printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
-           if ((dx <= dpx) && (dy <= dpy)) {
+           if (fDebugLevel>0)
+               printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
+           if ((dx <= dpx) && (dy <= dpy+eps)) {
                // consistent
                accepted[ico]=kTRUE;
                fprintf(stderr,"ico %d\n",ico);
@@ -700,8 +789,31 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 
        Float_t chi21 = 100;
        Float_t chi22 = 100;
+       Float_t chi23 = 100;
 
-       if (accepted[0]) {
+       fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
+           Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
+       
+       fQrInit[0]=fQrInit[1];
+
+       
+       if (accepted[0] && accepted[1]) {
+           fXInit[0]=xm[0][1];
+           fYInit[0]=0.5*(ym[0][0]+ym[0][1]);
+           fXInit[1]=xm[1][1];
+           fYInit[1]=0.5*(ym[0][0]+ym[1][1]);
+           fQrInit[0]=0.5;
+           fQrInit[1]=0.5;
+           chi23=CombiDoubleMathiesonFit(c);
+           if (chi23<10) {
+               Split(c);
+               Float_t yst;
+               yst = fYFit[0];
+               fYFit[0] = fYFit[1];
+               fYFit[1] = yst;
+               Split(c);
+           }
+       } else if (accepted[0]) {
            fXInit[0]=xm[0][0];
            fYInit[0]=ym[0][1];
            fXInit[1]=xm[1][1];
@@ -711,7 +823,8 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 //         Float_t prob = TMath::Prob(chi2,ndf);
 //         prob2->Fill(prob);
 //         chi2_2->Fill(chi21);
-           fprintf(stderr," chi2 %f\n",chi21);
+           if (fDebugLevel)
+               fprintf(stderr," chi2 %f\n",chi21);
            if (chi21<10) Split(c);
        } else if (accepted[1]) {
            fXInit[0]=xm[1][0];
@@ -723,11 +836,12 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 //         Float_t prob = TMath::Prob(chi2,ndf);
 //         prob2->Fill(prob);
 //         chi2_2->Fill(chi22);
-           fprintf(stderr," chi2 %f\n",chi22);
+           if (fDebugLevel)
+               fprintf(stderr," chi2 %f\n",chi22);
            if (chi22<10) Split(c);
        }
 
-       if (chi21 > 10 && chi22 > 10) {
+       if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
 //We keep only the combination found (X->cathode 2, Y->cathode 1)
            for (Int_t ico=0; ico<2; ico++) {
                if (accepted[ico]) {
@@ -736,6 +850,7 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
                    for (cath=0; cath<2; cath++) {
                        cnew.fX[cath]=Float_t(xm[ico][1]);
                        cnew.fY[cath]=Float_t(ym[ico][0]);
+                       cnew.fZ[cath]=fZPlane;
                        cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
                        for (i=0; i<fMul[cath]; i++) {
                            cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
@@ -756,7 +871,6 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
 //  (4) At least three local maxima on cathode 1 or on cathode 2 
 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     } else if (fNLocal[0]>2 || fNLocal[1]>2) {
-       
        Int_t param = fNLocal[0]*fNLocal[1];
        Int_t ii;
 
@@ -789,26 +903,30 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
        }
        
        Int_t nIco = ico;
-       
-       fprintf(stderr,"nIco %d\n",nIco);
+       if (fDebugLevel)
+           fprintf(stderr,"nIco %d\n",nIco);
        for (ico=0; ico<nIco; ico++) {
-           fprintf(stderr,"ico = %d\n",ico);
+           if (fDebugLevel)
+               fprintf(stderr,"ico = %d\n",ico);
            isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
            dpx=fSeg[0]->Dpx(isec)/2.;
            dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
            isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
            dpy=fSeg[1]->Dpy(isec)/2.;
            dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
-
-           fprintf(stderr,"dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy);
-           fprintf(stderr,"  X %f Y %f\n",xm[ico][1],ym[ico][0]);
+           if (fDebugLevel) {
+               fprintf(stderr,"dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy);
+               fprintf(stderr,"  X %f Y %f\n",xm[ico][1],ym[ico][0]);
+           }
            if ((dx <= dpx) && (dy <= dpy)) {
-               fprintf(stderr,"ok\n");
+               if (fDebugLevel)
+                   fprintf(stderr,"ok\n");
                Int_t cath;    
                AliMUONRawCluster cnew;
                for (cath=0; cath<2; cath++) {
                    cnew.fX[cath]=Float_t(xm[ico][1]);
                    cnew.fY[cath]=Float_t(ym[ico][0]);
+                   cnew.fZ[cath]=fZPlane;
                    cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
                    for (i=0; i<fMul[cath]; i++) {
                        cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
@@ -828,10 +946,11 @@ void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
     }
 }
 
-void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
+void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* /*c*/)
 {
 // Find all local maxima of a cluster
-    printf("\n Find Local maxima  !");
+    if (fDebugLevel)
+       printf("\n Find Local maxima  !");
     
     AliMUONDigit* digt;
     
@@ -867,12 +986,12 @@ void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
                digt=(AliMUONDigit*) fHitMap[cath]->GetHit(x[j], y[j]);
                isec=fSeg[cath]->Sector(x[j], y[j]);
                Float_t a1 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
-               if (digt->fSignal/a1 > fQ[i][cath]/a0) {
+               if (digt->Signal()/a1 > fQ[i][cath]/a0) {
                    isLocal[i][cath]=kFALSE;
                    break;
 //
 // handle special case of neighbouring pads with equal signal
-               } else if (digt->fSignal == fQ[i][cath]) {
+               } else if (digt->Signal() == fQ[i][cath]) {
                    if (fNLocal[cath]>0) {
                        for (Int_t k=0; k<fNLocal[cath]; k++) {
                            if (x[j]==fIx[fIndLocal[k][cath]][cath] 
@@ -890,11 +1009,13 @@ void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
            } 
        } // loop over all digits
     } // loop over cathodes
-    
-    printf("\n Found %d %d %d %d local Maxima\n",
-          fNLocal[0], fNLocal[1], fMul[0], fMul[1]);
-    fprintf(stderr,"\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]);
-    fprintf(stderr," Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]);
+
+    if (fDebugLevel) {
+       printf("\n Found %d %d %d %d local Maxima\n",
+              fNLocal[0], fNLocal[1], fMul[0], fMul[1]);
+       fprintf(stderr,"\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]);
+       fprintf(stderr," Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]);
+    }
     Int_t ix, iy, isec;
     Float_t dpx, dpy;
     
@@ -940,7 +1061,7 @@ void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
                if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
                    iNN++;
                    digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
-                   if (digt->fSignal > fQ[i][cath]) isLocal[i][cath]=kFALSE;
+                   if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
                }
            } // Loop over pad neighbours in y
            if (isLocal[i][cath] && iNN>0) {
@@ -950,12 +1071,14 @@ void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
        } // loop over all digits
 // if one additional maximum has been found we are happy 
 // if more maxima have been found restore the previous situation
-       fprintf(stderr,
-               "\n New search gives %d local maxima for cathode 1 \n",
-               fNLocal[0]);
-       fprintf(stderr,
-               "                  %d local maxima for cathode 2 \n",
-               fNLocal[1]);
+       if (fDebugLevel) {
+           fprintf(stderr,
+                   "\n New search gives %d local maxima for cathode 1 \n",
+                   fNLocal[0]);
+           fprintf(stderr,
+                   "                  %d local maxima for cathode 2 \n",
+                   fNLocal[1]);
+       }
        if (fNLocal[cath]>2) {
            fNLocal[cath]=iback;
        }
@@ -968,7 +1091,9 @@ void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
 //  Two local maxima on cathode 1 and one maximum on cathode 2 
 //  Look for local maxima considering left and right neighbours on the 2nd cathode only
        cath=1;
-       Int_t cath1=0;
+       Int_t cath1 = 0;
+       Float_t eps = 1.e-5;
+       
 //
 //  Loop over cluster digits
        for (i=0; i<fMul[cath]; i++) {
@@ -977,9 +1102,10 @@ void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
            dpy=fSeg[cath]->Dpy(isec);
            if (isLocal[i][cath]) continue;
 // Pad position should be consistent with position of local maxima on the opposite cathode
-           if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.) && 
-               (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.))
+           if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.+eps) && 
+               (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.+eps))
                continue;
+           
 //
 // get neighbours for that digit and assume that it is local maximum       
            isLocal[i][cath]=kTRUE;
@@ -988,22 +1114,23 @@ void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
 // iNN counts the number of neighbours with signal, it should be 1 or 2
            Int_t iNN=0;
            for (fSeg[cath]
-                    ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, 0., dpx);
+                    ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, dpx, 0.);
                 fSeg[cath]
                     ->MorePads();
                 fSeg[cath]
                     ->NextPad())
            {
+
                ix = fSeg[cath]->Ix();
                iy = fSeg[cath]->Iy();
-               
+
                // skip the current pad
                if (ix == fIx[i][cath]) continue;
                
                if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
                    iNN++;
                    digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
-                   if (digt->fSignal > fQ[i][cath]) isLocal[i][cath]=kFALSE;
+                   if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
                }
            } // Loop over pad neighbours in x
            if (isLocal[i][cath] && iNN>0) {
@@ -1013,9 +1140,11 @@ void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
        } // loop over all digits
 // if one additional maximum has been found we are happy 
 // if more maxima have been found restore the previous situation
-       fprintf(stderr,"\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]);
-       fprintf(stderr,"\n                  %d local maxima for cathode 2 \n",fNLocal[1]);
-//     printf("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]);
+       if (fDebugLevel) {
+           fprintf(stderr,"\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]);
+           fprintf(stderr,"\n                  %d local maxima for cathode 2 \n",fNLocal[1]);
+           printf("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]);
+       }
        if (fNLocal[cath]>2) {
            fNLocal[cath]=iback;
        }
@@ -1045,29 +1174,31 @@ void  AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_
        c->fQ[cath]=0;
     }
 
-//    fprintf(stderr,"\n fPeakSignal %d\n",c->fPeakSignal[cath]);
+    if (fDebugLevel)
+       fprintf(stderr,"\n fPeakSignal %d\n",c->fPeakSignal[cath]);
     for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
     {
        dig= fInput->Digit(cath,c->fIndexMap[i][cath]);
-       ix=dig->fPadX+c->fOffsetMap[i][cath];
-       iy=dig->fPadY;
-       Int_t q=dig->fSignal;
+       ix=dig->PadX()+c->fOffsetMap[i][cath];
+       iy=dig->PadY();
+       Int_t q=dig->Signal();
        if (!flag) q=Int_t(q*c->fContMap[i][cath]);
 //     fprintf(stderr,"q %d c->fPeakSignal[ %d ] %d\n",q,cath,c->fPeakSignal[cath]);
-       if (dig->fPhysics >= dig->fSignal) {
+       if (dig->Physics() >= dig->Signal()) {
            c->fPhysicsMap[i]=2;
-       } else if (dig->fPhysics == 0) {
+       } else if (dig->Physics() == 0) {
            c->fPhysicsMap[i]=0;
        } else  c->fPhysicsMap[i]=1;
 //
 // 
-//     fprintf(stderr,"q %d c->fPeakSignal[cath] %d\n",q,c->fPeakSignal[cath]);
+       if (fDebugLevel>1)
+           fprintf(stderr,"q %d c->fPeakSignal[cath] %d\n",q,c->fPeakSignal[cath]);
 // peak signal and track list
        if (q>c->fPeakSignal[cath]) {
            c->fPeakSignal[cath]=q;
-           c->fTracks[0]=dig->fHit;
-           c->fTracks[1]=dig->fTracks[0];
-           c->fTracks[2]=dig->fTracks[1];
+           c->fTracks[0]=dig->Hit();
+           c->fTracks[1]=dig->Track(0);
+           c->fTracks[2]=dig->Track(1);
 //         fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
        }
 //
@@ -1078,11 +1209,13 @@ void  AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_
            c->fQ[cath] += q;
        }
     } // loop over digits
-//    fprintf(stderr," fin du cluster c\n");
+    if (fDebugLevel)
+       fprintf(stderr," fin du cluster c\n");
 
 
     if (flag) {
        c->fX[cath]/=c->fQ[cath];
+// Force on anod
        c->fX[cath]=fSeg[cath]->GetAnod(c->fX[cath]);
        c->fY[cath]/=c->fQ[cath]; 
 //
@@ -1122,31 +1255,36 @@ void  AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath)
     {
        dig = fInput->Digit(cath,c->fIndexMap[i][cath]);
        fSeg[cath]->
-       GetPadC(dig->fPadX,dig->fPadY,xpad,ypad, zpad);
-       fprintf(stderr,"x %f y %f cx %f cy %f\n",xpad,ypad,c->fX[0],c->fY[0]);
+       GetPadC(dig->PadX(),dig->PadY(),xpad,ypad, zpad);
+       if (fDebugLevel)
+           fprintf(stderr,"x %f y %f cx %f cy %f\n",xpad,ypad,c->fX[0],c->fY[0]);
        dx = xpad - c->fX[0];
        dy = ypad - c->fY[0];
        dr = TMath::Sqrt(dx*dx+dy*dy);
 
        if (dr < dr0) {
            dr0 = dr;
-           fprintf(stderr," dr %f\n",dr);
-           Int_t q=dig->fSignal;
-           if (dig->fPhysics >= dig->fSignal) {
+           if (fDebugLevel)
+               fprintf(stderr," dr %f\n",dr);
+           Int_t q=dig->Signal();
+           if (dig->Physics() >= dig->Signal()) {
                c->fPhysicsMap[i]=2;
-           } else if (dig->fPhysics == 0) {
+           } else if (dig->Physics() == 0) {
                c->fPhysicsMap[i]=0;
            } else  c->fPhysicsMap[i]=1;
            c->fPeakSignal[cath]=q;
-           c->fTracks[0]=dig->fHit;
-           c->fTracks[1]=dig->fTracks[0];
-           c->fTracks[2]=dig->fTracks[1];
-           fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
+           c->fTracks[0]=dig->Hit();
+           c->fTracks[1]=dig->Track(0);
+           c->fTracks[2]=dig->Track(1);
+           if (fDebugLevel)
+               fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->Hit(),
+                   dig->Track(0));
        }
 //
     } // loop over digits
 
 //  apply correction to the coordinate along the anode wire
+// Force on anod
     c->fX[cath]=fSeg[cath]->GetAnod(c->fX[cath]);
 }
 
@@ -1162,15 +1300,15 @@ void  AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONR
     
     Int_t idx = fHitMap[cath]->GetHitIndex(i,j);
     AliMUONDigit* dig = (AliMUONDigit*) fHitMap[cath]->GetHit(i,j);
-    Int_t q=dig->fSignal;
-    Int_t theX=dig->fPadX;
-    Int_t theY=dig->fPadY
+    Int_t q=dig->Signal();
+    Int_t theX=dig->PadX();
+    Int_t theY=dig->PadY()
    
     if (q > TMath::Abs(c.fPeakSignal[0]) && q > TMath::Abs(c.fPeakSignal[1])) {
        c.fPeakSignal[cath]=q;
-       c.fTracks[0]=dig->fHit;
-       c.fTracks[1]=dig->fTracks[0];
-       c.fTracks[2]=dig->fTracks[1];
+       c.fTracks[0]=dig->Hit();
+       c.fTracks[1]=dig->Track(0);
+       c.fTracks[2]=dig->Track(1);
     }
 
 //
@@ -1179,9 +1317,9 @@ void  AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONR
     Int_t mu=c.fMultiplicity[cath];
     c.fIndexMap[mu][cath]=idx;
     
-    if (dig->fPhysics >= dig->fSignal) {
+    if (dig->Physics() >= dig->Signal()) {
         c.fPhysicsMap[mu]=2;
-    } else if (dig->fPhysics == 0) {
+    } else if (dig->Physics() == 0) {
         c.fPhysicsMap[mu]=0;
     } else  c.fPhysicsMap[mu]=1;
 
@@ -1189,9 +1327,9 @@ void  AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONR
     if (mu > 0) {
        for (Int_t ind = mu-1; ind >= 0; ind--) {
            Int_t ist=(c.fIndexMap)[ind][cath];
-           Int_t ql=fInput->Digit(cath, ist)->fSignal;
-           Int_t ix=fInput->Digit(cath, ist)->fPadX;
-           Int_t iy=fInput->Digit(cath, ist)->fPadY;
+           Int_t ql=fInput->Digit(cath, ist)->Signal();
+           Int_t ix=fInput->Digit(cath, ist)->PadX();
+           Int_t iy=fInput->Digit(cath, ist)->PadY();
            
            if (q>ql || (q==ql && theX > ix && theY < iy)) {
                c.fIndexMap[ind][cath]=idx;
@@ -1233,7 +1371,8 @@ void  AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONR
        iy=yList[in];
        
        if (fHitMap[cath]->TestHit(ix,iy)==kUnused) {
-//         printf("\n Neighbours %d %d %d", cath, ix, iy);
+           if (fDebugLevel>1)
+               printf("\n Neighbours %d %d %d", cath, ix, iy);
            FindCluster(ix, iy, cath, c);
        }
        
@@ -1263,11 +1402,13 @@ void  AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONR
     {
        
        ix = fSeg[iop]->Ix(); iy = fSeg[iop]->Iy();
-//         printf("\n ix, iy: %f %f %f %d %d %d", x,y,z,ix, iy, fSector);
+       if (fDebugLevel > 1)
+           printf("\n ix, iy: %f %f %f %d %d %d", x,y,z,ix, iy, fSector);
        if (fHitMap[iop]->TestHit(ix,iy)==kUnused){
            iXopp[nOpp]=ix;
            iYopp[nOpp++]=iy;
-//         printf("\n Opposite %d %d %d", iop, ix, iy);
+           if (fDebugLevel > 1)
+               printf("\n Opposite %d %d %d", iop, ix, iy);
        }
        
     } // Loop over pad neighbours
@@ -1311,45 +1452,50 @@ void AliMUONClusterFinderVS::FindRawClusters()
     for (cath=0; cath<2; cath++) {
        for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
            dig = fInput->Digit(cath, ndig);
-           Int_t i=dig->fPadX;
-           Int_t j=dig->fPadY;
+           Int_t i=dig->PadX();
+           Int_t j=dig->PadY();
            if (fHitMap[cath]->TestHit(i,j)==kUsed ||fHitMap[0]->TestHit(i,j)==kEmpty) {
                nskip++;
                continue;
            }
-           fprintf(stderr,"\n CATHODE %d CLUSTER %d\n",cath,ncls);
+           if (fDebugLevel)
+               fprintf(stderr,"\n CATHODE %d CLUSTER %d\n",cath,ncls);
            AliMUONRawCluster c;
            c.fMultiplicity[0]=0;
            c.fMultiplicity[1]=0;
-           c.fPeakSignal[cath]=dig->fSignal;
-           c.fTracks[0]=dig->fHit;
-           c.fTracks[1]=dig->fTracks[0];
-           c.fTracks[2]=dig->fTracks[1];
+           c.fPeakSignal[cath]=dig->Signal();
+           c.fTracks[0]=dig->Hit();
+           c.fTracks[1]=dig->Track(0);
+           c.fTracks[2]=dig->Track(1);
            // tag the beginning of cluster list in a raw cluster
            c.fNcluster[0]=-1;
            Float_t xcu, ycu;
            fSeg[cath]->GetPadC(i,j,xcu, ycu, fZPlane);
            fSector= fSeg[cath]->Sector(i,j)/100;
-//         printf("\n New Seed %d %d ", i,j);
+           if (fDebugLevel)
+               printf("\n New Seed %d %d ", i,j);
            
            FindCluster(i,j,cath,c);
 //          ^^^^^^^^^^^^^^^^^^^^^^^^
            // center of gravity
            c.fX[0] /= c.fQ[0];
+// Force on anod
            c.fX[0]=fSeg[0]->GetAnod(c.fX[0]);
            c.fY[0] /= c.fQ[0];
            c.fX[1] /= c.fQ[1];
+// Force on anod
            c.fX[1]=fSeg[0]->GetAnod(c.fX[1]);
            c.fY[1] /= c.fQ[1];
            
            c.fZ[0] = fZPlane;
            c.fZ[1] = fZPlane;      
 
-           fprintf(stderr,"\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
-                   c.fMultiplicity[0],c.fX[0],c.fY[0]);
-           fprintf(stderr," Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
-                   c.fMultiplicity[1],c.fX[1],c.fY[1]);
-//
+           if (fDebugLevel) {
+               fprintf(stderr,"\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
+                       c.fMultiplicity[0],c.fX[0],c.fY[0]);
+               fprintf(stderr," Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
+                       c.fMultiplicity[1],c.fX[1],c.fY[1]);
+           }
 //      Analyse cluster and decluster if necessary
 //     
        ncls++;
@@ -1413,7 +1559,9 @@ Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t c
     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
 // ready for minimisation      
-    clusterInput.Fitter()->SetPrintLevel(1);
+    clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
+    if (fDebugLevel==0)
+       clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
     clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
     arglist[0]= -1;
     arglist[1]= 0;
@@ -1440,7 +1588,7 @@ Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t c
     return fmin;
 }
 
-Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster *c)
+Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster * /*c*/)
 {
 // Perform combined Mathieson fit on both cathode planes
 //
@@ -1484,7 +1632,8 @@ Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster *c)
     if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
        
     icount=0;
-    printf("\n single y %f %f", fXInit[0], fYInit[0]);
+    if (fDebugLevel)
+       printf("\n single y %f %f", fXInit[0], fYInit[0]);
     
     for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy); 
         fSeg[0]->MorePads(); fSeg[0]->NextPad())
@@ -1493,7 +1642,8 @@ Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster *c)
        fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);     
        if (icount ==0) lower[1]=upper[1];
        icount++;
-       printf("\n upper lower %d %f %f", icount, upper[1], lower[1]);
+       if (fDebugLevel)
+           printf("\n upper lower %d %f %f", icount, upper[1], lower[1]);
     }
     
     if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
@@ -1504,7 +1654,9 @@ Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster *c)
     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
 // ready for minimisation      
-    clusterInput.Fitter()->SetPrintLevel(1);
+    clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
+    if (fDebugLevel==0)
+       clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
     clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
     arglist[0]= -1;
     arglist[1]= 0;
@@ -1531,7 +1683,7 @@ Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster *c)
     return fmin;
 }
 
-Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
+Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster * /*c*/, Int_t cath)
 {
 // Performs a double Mathieson fit on one cathode
 // 
@@ -1579,7 +1731,9 @@ Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster *c, Int_t ca
     clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
     clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
 // ready for minimisation      
-    clusterInput.Fitter()->SetPrintLevel(-1);
+    clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
+    if (fDebugLevel==0)
+       clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
     clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
     arglist[0]= -1;
     arglist[1]= 0;
@@ -1606,7 +1760,7 @@ Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster *c, Int_t ca
     return kTRUE;
 }
 
-Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster *c)
+Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster * /*c*/)
 {
 //
 // Perform combined double Mathieson fit on both cathode planes
@@ -1641,7 +1795,8 @@ Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster *c)
 
     Int_t icount;
     Float_t xdum, ydum, zdum;
-//    printf("\n Cluster Finder: %f %f %f %f  ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] );
+    if (fDebugLevel)
+       printf("\n Cluster Finder: %f %f %f %f  ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] );
     
 //  Find save upper and lower limits    
     icount = 0;
@@ -1650,22 +1805,30 @@ Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster *c)
         fSeg[1]->MorePads(); fSeg[1]->NextPad())
     {
        ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
+//     if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
        fSeg[1]->GetPadC(ix,iy,upper[0],ydum,zdum);     
        if (icount ==0) lower[0]=upper[0];
        icount++;
     }
     if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}    
+//    vstart[0] = 0.5*(lower[0]+upper[0]);
+
+    
     icount=0;
     
     for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy); 
         fSeg[0]->MorePads(); fSeg[0]->NextPad())
     {
        ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
+//     if (fHitMap[0]->TestHit(ix, iy) == kEmpty) continue;
        fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);     
        if (icount ==0) lower[1]=upper[1];
        icount++;
     }
+    
     if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}    
+//     vstart[1] = 0.5*(lower[1]+upper[1]);
+
 
     fSeg[1]->GetPadI(fXInit[1], fYInit[1], fZPlane, ix, iy);
     isec=fSeg[1]->Sector(ix, iy);
@@ -1683,11 +1846,13 @@ Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster *c)
         fSeg[1]->MorePads(); fSeg[1]->NextPad())
     {
        ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
+//     if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
        fSeg[1]->GetPadC(ix,iy,upper[2],ydum,zdum);     
        if (icount ==0) lower[2]=upper[2];
        icount++;
     }
     if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}    
+    //    vstart[2] = 0.5*(lower[2]+upper[2]);
 
     icount=0;
     
@@ -1695,12 +1860,17 @@ Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster *c)
         fSeg[0]-> MorePads(); fSeg[0]->NextPad())
     {
        ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
+//     if (fHitMap[0]->TestHit(ix, iy) != kEmpty) continue;
+       
        fSeg[0]->GetPadC(ix,iy,xdum,upper[3],zdum);     
        if (icount ==0) lower[3]=upper[3];
        icount++;
+
     }
     if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}    
-
+    
+//     vstart[3] = 0.5*(lower[3]+upper[3]);
+    
     lower[4]=0.;
     upper[4]=1.;
     lower[5]=0.;
@@ -1715,7 +1885,9 @@ Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster *c)
     clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
     clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
 // ready for minimisation      
-    clusterInput.Fitter()->SetPrintLevel(-1);
+    clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
+    if (fDebugLevel)
+       clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
     clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
     arglist[0]= -1;
     arglist[1]= 0;
@@ -1754,8 +1926,10 @@ void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
     for (j=0; j<2; j++) {
        AliMUONRawCluster cnew;
+       cnew.fGhost=c->fGhost;
        for (cath=0; cath<2; cath++) {
            cnew.fChi2[cath]=fChi2[0];
+           // ?? why not cnew.fChi2[cath]=fChi2[cath];
            
            if (fNPeaks == 0) {
                cnew.fNcluster[0]=-1;
@@ -1767,6 +1941,7 @@ void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
            cnew.fMultiplicity[cath]=0;
            cnew.fX[cath]=Float_t(fXFit[j]);
            cnew.fY[cath]=Float_t(fYFit[j]);
+           cnew.fZ[cath]=fZPlane;
            if (j==0) {
                cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]);
            } else {
@@ -1795,7 +1970,7 @@ void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
 //
 // Minimisation functions
 // Single Mathieson
-void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+void fcnS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
 {
     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
     Int_t i;
@@ -1815,7 +1990,7 @@ void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
     f=chisq;
 }
 
-void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+void fcnCombiS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
 {
     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
     Int_t i, cath;
@@ -1838,7 +2013,7 @@ void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t if
 }
 
 // Double Mathieson
-void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+void fcnS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
 {
     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
     Int_t i;
@@ -1860,7 +2035,7 @@ void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
 }
 
 // Double Mathieson
-void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+void fcnCombiS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
 {
     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
     Int_t i, cath;
@@ -1887,12 +2062,14 @@ void AliMUONClusterFinderVS::AddRawCluster(const AliMUONRawCluster c)
   // Add a raw cluster copy to the list
   //
     AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
-    pMUON->AddRawCluster(fInput->Chamber(),c); 
+    pMUON->GetMUONData()->AddRawCluster(fInput->Chamber(),c); 
     fNRawClusters++;
-    fprintf(stderr,"\nfNRawClusters %d\n",fNRawClusters);
+//    if (fDebugLevel)
+       fprintf(stderr,"\nfNRawClusters %d\n",fNRawClusters);
 }
 
 Bool_t AliMUONClusterFinderVS::TestTrack(Int_t t) {
+// Test if track was user selected
     if (fTrack[0]==-1 || fTrack[1]==-1) {
        return kTRUE;
     } else if (t==fTrack[0] || t==fTrack[1]) {
@@ -1903,7 +2080,7 @@ Bool_t AliMUONClusterFinderVS::TestTrack(Int_t t) {
 }
 
 AliMUONClusterFinderVS& AliMUONClusterFinderVS
-::operator = (const AliMUONClusterFinderVS& rhs)
+::operator = (const AliMUONClusterFinderVS& /*rhs*/)
 {
 // Dummy assignment operator
     return *this;