]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUON.cxx
Write a TreeD in SDigits2Digits method (needed to be compatible with alirun script)
[u/mrichter/AliRoot.git] / MUON / AliMUON.cxx
index 04a6585d5c9ab11c69ad035d549c2b6af7c3eab0..dbd928e756e295db4b4a8a9796d2f8acde131974 100644 (file)
  **************************************************************************/
 /*
 $Log$
+Revision 1.51  2001/05/31 10:19:52  morsch
+Fix for new AliRun::RunReco().
+
+Revision 1.50  2001/05/16 14:57:17  alibrary
+New files for folders and Stack
+
+Revision 1.49  2001/03/12 17:45:48  hristov
+Changes needed on Sun with CC 5.0
+
+Revision 1.48  2001/03/06 00:01:36  morsch
+Add  Digits2Reco() and FindClusters()
+Adapt call of cluster finder to new STEER.
+
+Revision 1.47  2001/03/05 08:38:36  morsch
+Digitization related methods moved to AliMUONMerger.
+
+Revision 1.46  2001/01/26 21:34:59  morsch
+Use access functions for AliMUONHit, AliMUONDigit and AliMUONPadHit data members.
+
+Revision 1.45  2001/01/26 20:00:49  hristov
+Major upgrade of AliRoot code
+
+Revision 1.44  2001/01/25 17:39:09  morsch
+Pass size of fNdch and fNrawch to CINT.
+
+Revision 1.43  2001/01/23 18:58:19  hristov
+Initialisation of some pointers
+
+Revision 1.42  2001/01/17 20:53:40  hristov
+Destructors corrected to avoid memory leaks
+
+Revision 1.41  2000/12/21 22:12:40  morsch
+Clean-up of coding rule violations,
+
+Revision 1.40  2000/11/29 20:32:26  gosset
+Digitize:
+1. correction for array index out of bounds
+2. one printout commented
+
+Revision 1.39  2000/11/12 17:17:03  pcrochet
+BuildGeometry of AliMUON for trigger chambers delegated to AliMUONSegmentationTriggerX (same strategy as for tracking chambers)
+
+Revision 1.38  2000/11/06 09:20:43  morsch
+AliMUON delegates part of BuildGeometry() to AliMUONSegmentation using the
+Draw() method. This avoids code and parameter replication.
+
+Revision 1.37  2000/10/26 09:53:37  pcrochet
+put back trigger chambers in the display (there was a problem in buildgeometry)
+
+Revision 1.36  2000/10/25 19:51:18  morsch
+Correct x-position of chambers.
+
+Revision 1.35  2000/10/24 19:46:21  morsch
+BuildGeometry updated for slats in station 3-4.
+
+Revision 1.34  2000/10/18 11:42:06  morsch
+- AliMUONRawCluster contains z-position.
+- Some clean-up of useless print statements during initialisations.
+
+Revision 1.33  2000/10/09 14:01:57  morsch
+Unused variables removed.
+
+Revision 1.32  2000/10/06 09:08:10  morsch
+Built geometry includes slat geometry for event display.
+
+Revision 1.31  2000/10/02 21:28:08  fca
+Removal of useless dependecies via forward declarations
+
+Revision 1.30  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.29  2000/07/28 13:49:38  morsch
+SetAcceptance defines inner and outer chamber radii according to angular acceptance.
+Can be used for simple acceptance studies.
+
+Revision 1.28  2000/07/22 16:43:15  morsch
+Same comment as before, but now done correctly I hope (sorry it's Saturday evening)
+
+Revision 1.27  2000/07/22 16:36:50  morsch
+Change order of indices in creation (new) of xhit and yhit
+
+Revision 1.26  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.25  2000/06/29 12:34:09  morsch
+AliMUONSegmentation class has been made independent of AliMUONChamber. This makes
+it usable with any other geometry class. The link to the object to which it belongs is
+established via an index. This assumes that there exists a global geometry manager
+from which the pointer to the parent object can be obtained (in our case gAlice).
+
+Revision 1.24  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.23  2000/06/28 12:19:17  morsch
+More consequent seperation of global input data services (AliMUONClusterInput singleton) and the
+cluster and hit reconstruction algorithms in AliMUONClusterFindRawinderVS.
+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.22  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.
@@ -97,6 +206,7 @@ Log message added
 #include <TTUBE.h>
 #include <TBRIK.h>
 #include <TRotMatrix.h>
+#include <TGeometry.h>
 #include <TNode.h> 
 #include <TTree.h> 
 #include <TRandom.h> 
@@ -123,15 +233,18 @@ Log message added
 #include "AliMUONRawCluster.h"
 #include "AliMUONLocalTrigger.h"
 #include "AliMUONGlobalTrigger.h"
-#include "AliMUONHitMap.h"
+#include "AliMUONTriggerCircuit.h"
+#include "AliHitMap.h"
 #include "AliMUONHitMapA1.h"
 #include "AliMUONChamberTrigger.h"
 #include "AliMUONConstants.h"
 #include "AliMUONClusterFinderVS.h"
 #include "AliMUONTriggerDecision.h"
 #include "AliRun.h"
+#include "AliHeader.h"
 #include "AliMC.h"
 #include "AliMUONClusterInput.h"
+#include "AliMUONMerger.h"     
 #include "iostream.h"
 #include "AliCallf77.h" 
 #include "AliConst.h" 
@@ -154,19 +267,28 @@ ClassImp(AliMUON)
 //___________________________________________
 AliMUON::AliMUON()
 {
-   fIshunt       = 0;
-   fHits         = 0;
-   fPadHits      = 0;
-   fNPadHits     = 0;
-   fDchambers    = 0;
-   fTriggerCircuits = 0;     // cp new design of AliMUONTriggerDecision
-   fNdch         = 0;
-   fRawClusters  = 0;
-   fNrawch       = 0;
-   fGlobalTrigger   = 0;
-   fNLocalTrigger   = 0;
-   fLocalTrigger    = 0;
-   fNLocalTrigger   = 0;
+// Default Constructor
+//
+    fNCh             = 0;
+    fNTrackingCh     = 0;
+    fIshunt          = 0;
+    fHits            = 0;
+    fPadHits         = 0;
+    fNPadHits        = 0;
+    fChambers        = 0;
+    fDchambers       = 0;
+    fTriggerCircuits = 0;  
+    fNdch            = 0;
+    fRawClusters     = 0;
+    fNrawch          = 0;
+    fGlobalTrigger   = 0;
+    fNLocalTrigger   = 0;
+    fLocalTrigger    = 0;
+    fNLocalTrigger   = 0;
+    fAccMin          = 0.;
+    fAccMax          = 0.;   
+    fAccCut          = kFALSE;
+    fMerger          = 0;
 }
  
 //___________________________________________
@@ -185,7 +307,9 @@ AliMUON::AliMUON(const char *name, const char *title)
    fNPadHits  =  0;
    fIshunt     =  0;
 
-   fNdch      = new Int_t[AliMUONConstants::NCh()];
+   fNCh             = AliMUONConstants::NCh(); 
+   fNTrackingCh     = AliMUONConstants::NTrackingCh();
+   fNdch            = new Int_t[fNCh];
 
    fDchambers = new TObjArray(AliMUONConstants::NCh());
 
@@ -196,7 +320,7 @@ AliMUON::AliMUON(const char *name, const char *title)
        fNdch[i]=0;
    }
 
-   fNrawch      = new Int_t[AliMUONConstants::NTrackingCh()];
+   fNrawch      = new Int_t[fNTrackingCh];
 
    fRawClusters = new TObjArray(AliMUONConstants::NTrackingCh());
 
@@ -204,19 +328,12 @@ AliMUON::AliMUON(const char *name, const char *title)
        (*fRawClusters)[i] = new TClonesArray("AliMUONRawCluster",10000); 
        fNrawch[i]=0;
    }
-   cout << " here " << "\n";
 
    fGlobalTrigger = new TClonesArray("AliMUONGlobalTrigger",1);    
    fNGlobalTrigger = 0;
    fLocalTrigger  = new TClonesArray("AliMUONLocalTrigger",234);    
    fNLocalTrigger = 0;
 
-//   
-// Transport angular cut
-   fAccCut=0;
-   fAccMin=2;
-   fAccMax=9;
-
    SetMarkerColor(kRed);
 //
 //
@@ -239,9 +356,9 @@ AliMUON::AliMUON(const char *name, const char *title)
            ch = 2 * st + stCH;
 //
            if (ch < AliMUONConstants::NTrackingCh()) {
-               (*fChambers)[ch] = new AliMUONChamber();
+               (*fChambers)[ch] = new AliMUONChamber(ch);
            } else {
-               (*fChambers)[ch] = new AliMUONChamberTrigger();
+               (*fChambers)[ch] = new AliMUONChamberTrigger(ch);
            }
            
            AliMUONChamber* chamber = (AliMUONChamber*) (*fChambers)[ch];
@@ -251,6 +368,7 @@ AliMUON::AliMUON(const char *name, const char *title)
            chamber->SetZ(AliMUONConstants::DefaultChamberZ(ch));
 //
            chamber->InitGeo(AliMUONConstants::DefaultChamberZ(ch));
+//          Set chamber inner and outer radius to default
            chamber->SetRInner(AliMUONConstants::Dmin(st)/2);
            chamber->SetROuter(AliMUONConstants::Dmax(st)/2);
 //
@@ -263,14 +381,18 @@ AliMUON::AliMUON(const char *name, const char *title)
 //
    fMaxIterPad   = 0;
    fCurIterPad   = 0;
+//
+   fAccMin          = 0.;
+   fAccMax          = 0.;   
+   fAccCut          = kFALSE;
 
    // cp new design of AliMUONTriggerDecision
    fTriggerCircuits = new TObjArray(AliMUONConstants::NTriggerCircuit());
    for (Int_t circ=0; circ<AliMUONConstants::NTriggerCircuit(); circ++) {
      (*fTriggerCircuits)[circ] = new AliMUONTriggerCircuit();     
-   }
-   // cp new design of AliMUONTriggerDecision
 
+   }
+     fMerger = 0;
 }
  
 //___________________________________________
@@ -283,35 +405,79 @@ AliMUON::AliMUON(const AliMUON& rMUON)
 
 AliMUON::~AliMUON()
 {
-    printf("Calling AliMUON destructor !!!\n");
+// Destructor
+    if(fDebug) printf("%s: Calling AliMUON destructor !!!\n",ClassName());
     
     Int_t i;
     fIshunt  = 0;
-    delete fHits;
-    delete fPadHits;
-    
-    delete fGlobalTrigger;
+    // Delete TObjArrays
+    if (fChambers){
+      fChambers->Delete();
+      delete fChambers;
+    }
+    if (fTriggerCircuits){
+      fTriggerCircuits->Delete();
+      delete fTriggerCircuits;
+    }
+    if (fDchambers){
+      fDchambers->Delete();
+      delete fDchambers;
+    }
+    if (fRawClusters){
+      fRawClusters->Delete();
+      delete fRawClusters;
+    }
+    for (i=0;i<AliMUONConstants::NTrackingCh();i++) {
+      fNrawch[i]=0;
+    }
+    // Delete TClonesArrays
+    if (fPadHits){
+      fPadHits->Delete();
+      delete fPadHits;
+    }
+    if (fGlobalTrigger){
+      fGlobalTrigger->Delete();
+      delete fGlobalTrigger;
+    }
     fNGlobalTrigger = 0;
     
-    delete fLocalTrigger;
+    if (fLocalTrigger){
+      fLocalTrigger->Delete();
+      delete fLocalTrigger;
+    }
     fNLocalTrigger = 0;
 
-    for (i=0;i<AliMUONConstants::NCh();i++) {
-       delete (*fDchambers)[i];
-       fNdch[i]=0;
+    if (fHits2){
+      fHits2->Delete();
+      delete fHits2;
     }
-    delete fDchambers;
-    
-    for (i=0;i<AliMUONConstants::NTrackingCh();i++) {
-       delete (*fRawClusters)[i];
-       fNrawch[i]=0;
+
+    if (fPadHits2){
+      fPadHits2->Delete();
+      delete fPadHits2;
     }
-    delete fRawClusters;
-    
-    for (Int_t circ=0; circ<AliMUONConstants::NTriggerCircuit(); circ++) {
-       delete (*fTriggerCircuits)[circ];
+
+    if (fHits) {
+      fHits->Delete();
+      delete fHits;
+    }
+
+    // Delete hits tree for background event
+
+    if (fTrH1) {
+      fTrH1->Delete();
+      delete fTrH1;
     }
-    delete fTriggerCircuits;
+
+    if (fMerger) delete fMerger;
 }
  
 //___________________________________________
@@ -370,175 +536,15 @@ void AliMUON::AddLocalTrigger(Int_t *localtr)
 //___________________________________________
 void AliMUON::BuildGeometry()
 {
-    TNode *node, *nodeF, *top, *nodeS;
-    const int kColorMUON  = kBlue;
-    const int kColorMUON2 = kGreen; 
-    //
-    top=gAlice->GetGeometry()->GetNode("alice");
-// MUON
-//
-//  z-Positions of Chambers
-    const Float_t kCz[7]={511., 686., 971., 1245., 1445., 1600, 1700.};
-//  inner diameter (Xlenght for trigger chamber -> active area)
-    const Float_t kDmin[7]={ 35.,  47.,  67.,   86.,  100., 544., 544.};
-//  outer diameter (Ylenght for trigger chamber -> active area)
-    const Float_t kDmax[7]={183., 245., 346.,  520.,  520., 612., 612.};
-
-    TRotMatrix* rot000 = new TRotMatrix("Rot000"," ", 90,  0, 90, 90, 0, 0);
-    TRotMatrix* rot090 = new TRotMatrix("Rot090"," ", 90, 90, 90,180, 0, 0);
-    TRotMatrix* rot180 = new TRotMatrix("Rot180"," ", 90,180, 90,270, 0, 0);
-    TRotMatrix* rot270 = new TRotMatrix("Rot270"," ", 90,270, 90,  0, 0, 0);
-    
-    Float_t rmin, rmax, dx, dy, dz, dr, xpos, ypos, zpos;
-    Float_t dzc1=4.;           // tracking chambers
-    Float_t dzc2=15.;          // trigger chambers
-    Float_t hole=102.;          // x-y hole around beam pipe for trig. chambers
-    Float_t zscale;            // scaling parameter trigger chambers
-    Float_t halfx, halfy;   
-    char nameChamber[9], nameSense[9], nameFrame[9], nameNode[8];
-    char nameSense1[9], nameSense2[9];    
-    for (Int_t i=0; i<7; i++) {
-       for (Int_t j=0; j<2; j++) {
-           Int_t id=2*i+j+1;
-           if (i<5) {               // tracking chambers
-               if (j==0) {
-                   zpos=kCz[i]-dzc1;
-               } else {
-                   zpos=kCz[i]+dzc1;
-               }
-           } else {
-               if (j==0) {
-                   zpos=kCz[i];
-               } else {            
-                   zpos=kCz[i]+dzc2;
-               }
-           }
-           sprintf(nameChamber,"C_MUON%d",id);
-           sprintf(nameSense,"S_MUON%d",id);
-           sprintf(nameSense1,"S1_MUON%d",id);
-           sprintf(nameSense2,"S2_MUON%d",id);
-           sprintf(nameFrame,"F_MUON%d",id);   
-           if (i<5) {                        // tracking chambers
-               rmin = kDmin[i]/2.-3;
-               rmax = kDmax[i]/2.+3;
-               new TTUBE(nameChamber,"Mother","void",rmin,rmax,0.25,1.);
-               rmin = kDmin[i]/2.;
-               rmax = kDmax[i]/2.;
-               new TTUBE(nameSense,"Sens. region","void",rmin,rmax,0.25, 1.);
-               dx=(rmax-rmin)/2;
-               dy=3.;
-               dz=0.25;
-               TBRIK* frMUON = new TBRIK(nameFrame,"Frame","void",dx,dy,dz);
-               top->cd();
-               sprintf(nameNode,"MUON%d",100+id);
-               node = new TNode(nameNode,"ChamberNode",nameChamber,0,0,zpos,"");
-               node->SetLineColor(kColorMUON);
-               fNodes->Add(node);
-               node->cd();
-               sprintf(nameNode,"MUON%d",200+id);
-               node = new TNode(nameNode,"Sens. Region Node",nameSense,0,0,0,"");
-               node->SetLineColor(kColorMUON);
-               node->cd();
-               dr=dx+rmin;
-               sprintf(nameNode,"MUON%d",300+id);
-               nodeF = new TNode(nameNode,"Frame0",frMUON,dr, 0, 0,rot000,"");
-               nodeF->SetLineColor(kColorMUON);
-               node->cd();
-               sprintf(nameNode,"MUON%d",400+id);
-               nodeF = new TNode(nameNode,"Frame1",frMUON,0 ,dr,0,rot090,"");
-               nodeF->SetLineColor(kColorMUON);
-               node->cd();
-               sprintf(nameNode,"MUON%d",500+id);
-               nodeF = new TNode(nameNode,"Frame2",frMUON,-dr,0,0,rot180,"");
-               nodeF->SetLineColor(kColorMUON);
-               node  ->cd();
-               sprintf(nameNode,"MUON%d",600+id);   
-               nodeF = new TNode(nameNode,"Frame3",frMUON,0,-dr,0,rot270,"");
-               nodeF->SetLineColor(kColorMUON);
-           } else { 
-               zscale=zpos/kCz[5];
-               Float_t xsize=kDmin[i]*zscale;
-               Float_t ysize=kDmax[i]*zscale;
-               Float_t holeScaled=hole*zscale;
-               
-               halfx=xsize/2.+3.;
-               halfy=ysize/2.+3.;          
-               new TBRIK(nameChamber,"Mother","void",halfx,halfy,0.25);
-               top->cd();
-               sprintf(nameNode,"MUON%d",100+id);
-               node = new TNode(nameNode,"Chambernode",nameChamber,0,0,zpos,"");
-               node->SetLineColor(kColorMUON2);
-               fNodes->Add(node);
-               
-// up/down of beam pipe
-               halfx=xsize/2.;
-               halfy=(ysize/2.-holeScaled/2.)/2.;          
-               new TBRIK(nameSense,"Sens. region","void",halfx,halfy,0.25);
-               
-               node->cd();
-               ypos=holeScaled/2.+((ysize/2.-holeScaled/2.)/2.);
-               sprintf(nameNode,"MUON%d",200+id);
-               nodeS = new TNode(nameNode,"Sens. Region Node",nameSense,0,ypos,0,"");
-               nodeS->SetLineColor(kColorMUON2);
-               
-               node->cd();
-               ypos=-1.*ypos;
-               sprintf(nameNode,"MUON%d",300+id);
-               nodeS = new TNode(nameNode,"Sens. Region Node",nameSense,0,ypos,0,"");
-               nodeS->SetLineColor(kColorMUON2);
-               
-// left/right of beam pipe
-               halfx=(xsize/2.-holeScaled/2.)/2.;
-               halfy=holeScaled/2.;    
-               new TBRIK(nameSense1,"Sens. region","void",halfx,halfy,0.25);
-               
-               node->cd();
-               xpos=holeScaled/2.+((xsize/2.-holeScaled/2.)/2.);           
-               sprintf(nameNode,"MUON%d",400+id);
-               nodeS = new TNode(nameNode,"Sens. Region Node",nameSense1,xpos,0,0,"");
-               nodeS->SetLineColor(kColorMUON2);
-               
-               node->cd();
-               xpos=-1.*xpos;
-               sprintf(nameNode,"MUON%d",500+id);
-               nodeS = new TNode(nameNode,"Sens. Region Node",nameSense1,xpos,0,0,"");
-               nodeS->SetLineColor(kColorMUON2);
-               
-// missing corners
-               halfx=17.*zscale/2.;
-               halfy=halfx;
-               new TBRIK(nameSense2,"Sens. region","void",halfx,halfy,0.25);
-               
-               node->cd();
-               xpos=holeScaled/2.-halfx;
-               ypos=xpos;
-               sprintf(nameNode,"MUON%d",600+id);
-               nodeS = new TNode(nameNode,"Sens. Region Node",nameSense2,xpos,ypos,0,"");
-               nodeS->SetLineColor(kColorMUON2);
-               
-               node->cd();
-               ypos=-1.*xpos;
-               sprintf(nameNode,"MUON%d",700+id);
-               nodeS = new TNode(nameNode,"Sens. Region Node",nameSense2,xpos,ypos,0,"");
-               nodeS->SetLineColor(kColorMUON2);
-               
-               node->cd();
-               xpos=-1.*xpos;
-               sprintf(nameNode,"MUON%d",800+id);
-               nodeS = new TNode(nameNode,"Sens. Region Node",nameSense2,xpos,ypos,0,"");
-               nodeS->SetLineColor(kColorMUON2);
-               
-               node->cd();
-               ypos=-1.*xpos;
-               sprintf(nameNode,"MUON%d",900+id);
-               nodeS = new TNode(nameNode,"Sens. Region Node",nameSense2,xpos,ypos,0,"");
-               nodeS->SetLineColor(kColorMUON2);
-           } 
-       }
+// Geometry for event display
+  for (Int_t i=0; i<7; i++) {
+    for (Int_t j=0; j<2; j++) {
+      Int_t id=2*i+j+1;
+      this->Chamber(id-1).SegmentationModel(1)->Draw("eventdisplay");
     }
+  }
 }
 
-
 //___________________________________________
 Int_t AliMUON::DistancetoPrimitive(Int_t , Int_t )
 {
@@ -546,57 +552,77 @@ Int_t AliMUON::DistancetoPrimitive(Int_t , Int_t )
 }
 
 //___________________________________________
-void AliMUON::MakeBranch(Option_t* option)
+void AliMUON::MakeBranch(Option_t* option, const char *file)
 {
+    //
     // Create Tree branches for the MUON.
+    //
     const Int_t kBufferSize = 4000;
     char branchname[30];
     sprintf(branchname,"%sCluster",GetName());
     
-    AliDetector::MakeBranch(option);
-    
-    if (fPadHits   && gAlice->TreeH()) {
-       gAlice->TreeH()->Branch(branchname,&fPadHits, kBufferSize);
-       printf("Making Branch %s for clusters\n",branchname);
-    }
-    
-// one branch for digits per chamber
-    Int_t i;
+    AliDetector::MakeBranch(option,file);
     
-    for (i=0; i<AliMUONConstants::NCh() ;i++) {
-       sprintf(branchname,"%sDigits%d",GetName(),i+1);
-       
-       if (fDchambers   && gAlice->TreeD()) {
-           gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), kBufferSize);
-           printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
-       }       
+    const char *cD = strstr(option,"D");
+    const char *cR = strstr(option,"R");
+    const char *cH = strstr(option,"H");
+
+    if (fPadHits   && gAlice->TreeH() && cH) {
+      MakeBranchInTree(gAlice->TreeH(), 
+                       branchname, &fPadHits, kBufferSize, file);
     }
     
-    printf("Make Branch - TreeR address %p\n",gAlice->TreeR());
+    if (cD) {
+      //
+      // one branch for digits per chamber
+      // 
+      Int_t i;
     
-// one branch for raw clusters per chamber
-    for (i=0; i<AliMUONConstants::NTrackingCh() ;i++) {
-       sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
-       
-       if (fRawClusters   && gAlice->TreeR()) {
-           gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), kBufferSize);
-           printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
-       }       
-    }
-    
-// one branch for global trigger
-    sprintf(branchname,"%sGlobalTrigger",GetName());
-    if (fGlobalTrigger && gAlice->TreeR()) {  
-       gAlice->TreeR()->Branch(branchname,&fGlobalTrigger,kBufferSize);
-       printf("Making Branch %s for Global Trigger\n",branchname);
-    }
-// one branch for local trigger
-    sprintf(branchname,"%sLocalTrigger",GetName());
-    if (fLocalTrigger && gAlice->TreeR()) {  
-       gAlice->TreeR()->Branch(branchname,&fLocalTrigger,kBufferSize);
-       printf("Making Branch %s for Local Trigger\n",branchname);
+      for (i=0; i<AliMUONConstants::NCh() ;i++) {
+           sprintf(branchname,"%sDigits%d",GetName(),i+1);     
+           if (fDchambers   && gAlice->TreeD()) {
+            MakeBranchInTree(gAlice->TreeD(), 
+                             branchname, &((*fDchambers)[i]), kBufferSize, file);
+             printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
+        }
+         }     
     }
     
+    if (cR) {
+      //     
+      // one branch for raw clusters per chamber
+      //  
+      printf("Make Branch - TreeR address %p\n",gAlice->TreeR());
+      
+      Int_t i;
+
+      for (i=0; i<AliMUONConstants::NTrackingCh() ;i++) {
+           sprintf(branchname,"%sRawClusters%d",GetName(),i+1);        
+           if (fRawClusters   && gAlice->TreeR()) {
+              MakeBranchInTree(gAlice->TreeR(), 
+                               branchname, &((*fRawClusters)[i]), kBufferSize, file);
+             printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
+           }   
+      }
+      //
+      // one branch for global trigger
+      //
+      sprintf(branchname,"%sGlobalTrigger",GetName());
+      if (fGlobalTrigger && gAlice->TreeR()) {  
+          MakeBranchInTree(gAlice->TreeR(), 
+                           branchname, &fGlobalTrigger, kBufferSize, file);
+           printf("Making Branch %s for Global Trigger\n",branchname);
+      }
+      //
+      // one branch for local trigger
+      //  
+      sprintf(branchname,"%sLocalTrigger",GetName());
+      if (fLocalTrigger && gAlice->TreeR()) {  
+          MakeBranchInTree(gAlice->TreeR(), 
+                           branchname, &fLocalTrigger, kBufferSize, file);
+           printf("Making Branch %s for Local Trigger\n",branchname);
+      }
+   }
 }
 
 //___________________________________________
@@ -694,6 +720,7 @@ void AliMUON::ResetTrigger()
 //____________________________________________
 void AliMUON::SetPadSize(Int_t id, Int_t isec, Float_t p1, Float_t p2)
 {
+// Set the pad size for chamber id and cathode isec
     Int_t i=2*(id-1);
     ((AliMUONChamber*) (*fChambers)[i])  ->SetPadSize(isec,p1,p2);
     ((AliMUONChamber*) (*fChambers)[i+1])->SetPadSize(isec,p1,p2);
@@ -721,6 +748,7 @@ void AliMUON::SetChambersZToDefault()
 //___________________________________________
 void AliMUON::SetChargeSlope(Int_t id, Float_t p1)
 {
+// Set the inverse charge slope for chamber id
     Int_t i=2*(id-1);
     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSlope(p1);
     ((AliMUONChamber*) (*fChambers)[i+1])->SetChargeSlope(p1);
@@ -729,6 +757,7 @@ void AliMUON::SetChargeSlope(Int_t id, Float_t p1)
 //___________________________________________
 void AliMUON::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
 {
+// Set sigma of charge spread for chamber id
     Int_t i=2*(id-1);
     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSpread(p1,p2);
     ((AliMUONChamber*) (*fChambers)[i+1])->SetChargeSpread(p1,p2);
@@ -737,6 +766,7 @@ void AliMUON::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
 //___________________________________________
 void AliMUON::SetSigmaIntegration(Int_t id, Float_t p1)
 {
+// Set integration limits for charge spread
     Int_t i=2*(id-1);
     ((AliMUONChamber*) (*fChambers)[i])->SetSigmaIntegration(p1);
     ((AliMUONChamber*) (*fChambers)[i+1])->SetSigmaIntegration(p1);
@@ -745,6 +775,7 @@ void AliMUON::SetSigmaIntegration(Int_t id, Float_t p1)
 //___________________________________________
 void AliMUON::SetMaxAdc(Int_t id, Int_t p1)
 {
+// Set maximum number for ADCcounts (saturation)
     Int_t i=2*(id-1);
     ((AliMUONChamber*) (*fChambers)[i])->SetMaxAdc(p1);
     ((AliMUONChamber*) (*fChambers)[i+1])->SetMaxAdc(p1);
@@ -753,61 +784,99 @@ void AliMUON::SetMaxAdc(Int_t id, Int_t p1)
 //___________________________________________
 void AliMUON::SetMaxStepGas(Float_t p1)
 {
+// Set stepsize in gas
      fMaxStepGas=p1;
 }
 
 //___________________________________________
 void AliMUON::SetMaxStepAlu(Float_t p1)
 {
+// Set step size in Alu
     fMaxStepAlu=p1;
 }
 
 //___________________________________________
 void AliMUON::SetMaxDestepGas(Float_t p1)
 {
+// Set maximum step size in Gas
     fMaxDestepGas=p1;
 }
 
 //___________________________________________
 void AliMUON::SetMaxDestepAlu(Float_t p1)
 {
+// Set maximum step size in Alu
     fMaxDestepAlu=p1;
 }
 //___________________________________________
-void AliMUON::SetMuonAcc(Bool_t acc, Float_t angmin, Float_t angmax)
+void AliMUON::SetAcceptance(Bool_t acc, Float_t angmin, Float_t angmax)
 {
+// Set acceptance cuts 
    fAccCut=acc;
-   fAccMin=angmin;
-   fAccMax=angmax;
+   fAccMin=angmin*TMath::Pi()/180;
+   fAccMax=angmax*TMath::Pi()/180;
+   Int_t ch;
+   if (acc) {
+       for (Int_t st = 0; st < AliMUONConstants::NCh() / 2; st++) {
+          // Loop over 2 chambers in the station
+          for (Int_t stCH = 0; stCH < 2; stCH++) {
+              ch = 2 * st + stCH;
+//         Set chamber inner and outer radius according to acceptance cuts
+              Chamber(ch).SetRInner(AliMUONConstants::DefaultChamberZ(ch)*TMath::Tan(fAccMin));
+              Chamber(ch).SetROuter(AliMUONConstants::DefaultChamberZ(ch)*TMath::Tan(fAccMax));
+          } // chamber loop
+       } // station loop
+   }
 }
 //___________________________________________
-void   AliMUON::SetSegmentationModel(Int_t id, Int_t isec, AliMUONSegmentation *segmentation)
+void   AliMUON::SetSegmentationModel(Int_t id, Int_t isec, AliSegmentation *segmentation)
 {
+// Set the segmentation for chamber id cathode isec
     ((AliMUONChamber*) (*fChambers)[id])->SetSegmentationModel(isec, segmentation);
 
 }
 //___________________________________________
 void   AliMUON::SetResponseModel(Int_t id, AliMUONResponse *response)
 {
+// Set the response for chamber id
     ((AliMUONChamber*) (*fChambers)[id])->SetResponseModel(response);
 }
 
 void   AliMUON::SetReconstructionModel(Int_t id, AliMUONClusterFinderVS *reconst)
 {
+// Set ClusterFinder for chamber id
     ((AliMUONChamber*) (*fChambers)[id])->SetReconstructionModel(reconst);
 }
 
 void   AliMUON::SetNsec(Int_t id, Int_t nsec)
 {
+// Set number of segmented cathods for chamber id
     ((AliMUONChamber*) (*fChambers)[id])->SetNsec(nsec);
 }
 
-
 //___________________________________________
+void AliMUON::SDigits2Digits()
+{
 
+// write TreeD here 
 
+    if (!fMerger) {
+      if (gAlice->GetDebug()>0) {
+       cerr<<"AliMUON::SDigits2Digits: create default AliMUONMerger "<<endl;
+       cerr<<" no merging, just digitization of 1 event will be done"<<endl;
+      }
+      fMerger = new AliMUONMerger();
+    }
+    fMerger->Init();
+    fMerger->Digitise();
+    char hname[30];
+    sprintf(hname,"TreeD%d",gAlice->GetHeader()->GetEvent());
+    gAlice->TreeD()->Write(hname,TObject::kOverwrite);
+    gAlice->TreeD()->Reset();
+}
 
-void AliMUON::MakePadHits(Float_t xhit,Float_t yhit,
+//___________________________________________
+void AliMUON::MakePadHits(Float_t xhit,Float_t yhit, Float_t zhit,
                          Float_t eloss, Float_t tof,  Int_t idvol)
 {
 //
@@ -826,10 +895,13 @@ void AliMUON::MakePadHits(Float_t xhit,Float_t yhit,
     clhits[0]=fNhits+1;
 //
 //
+//    if (idvol == 6) printf("\n ->Disintegration %f %f %f", xhit, yhit, eloss );
+    
+
     ((AliMUONChamber*) (*fChambers)[idvol])
-       ->DisIntegration(eloss, tof, xhit, yhit, nnew, newclust);
+       ->DisIntegration(eloss, tof, xhit, yhit, zhit, nnew, newclust);
     Int_t ic=0;
-    
+//    if (idvol == 6) printf("\n nnew  %d \n", nnew);
 //
 //  Add new clusters
     for (Int_t i=0; i<nnew; i++) {
@@ -853,482 +925,6 @@ void AliMUON::MakePadHits(Float_t xhit,Float_t yhit,
     }
 }
 
-//----------------------------------------------------------------------
-
-void AliMUON::Digitise(Int_t nev,Int_t bgrEvent,Option_t *option,Option_t *opt,Text_t *filename)
-{
-    // keep galice.root for signal and name differently the file for 
-    // background when add! otherwise the track info for signal will be lost !
-  
-    static Bool_t first=kTRUE;
-    static TFile *file;
-    char *addBackground = strstr(option,"Add");
-
-
-    AliMUONChamber*  iChamber;
-    AliMUONSegmentation*  segmentation;
-
-    
-    Int_t trk[50];
-    Int_t chtrk[50];  
-    TObjArray *list=new TObjArray;
-    static TClonesArray *pAddress=0;
-    if(!pAddress) pAddress=new TClonesArray("TVector",1000);
-    Int_t digits[5]; 
-
-    AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
-    AliMUONHitMap** hitMap= new AliMUONHitMap* [AliMUONConstants::NCh()];
-    for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {hitMap[i]=0;}
-    if (addBackground ) {
-       if(first) {
-           fFileName=filename;
-           cout<<"filename"<<fFileName<<endl;
-           file=new TFile(fFileName);
-           cout<<"I have opened "<<fFileName<<" file "<<endl;
-           fHits2     = new TClonesArray("AliMUONHit",1000  );
-           fPadHits2 = new TClonesArray("AliMUONPadHit",10000);
-       }           
-       first=kFALSE;
-       file->cd();
-       //file->ls();
-       // Get Hits Tree header from file
-       if(fHits2) fHits2->Clear();
-       if(fPadHits2) fPadHits2->Clear();
-       if(fTrH1) delete fTrH1;
-       fTrH1=0;
-       
-       char treeName[20];
-       sprintf(treeName,"TreeH%d",bgrEvent);
-       fTrH1 = (TTree*)gDirectory->Get(treeName);
-        //printf("fTrH1 %p of treename %s for event %d \n",fTrH1,treeName,bgrEvent);
-       
-       if (!fTrH1) {
-           printf("ERROR: cannot find Hits Tree for event:%d\n",bgrEvent);
-       }
-       // Set branch addresses
-       TBranch *branch;
-       char branchname[20];
-       sprintf(branchname,"%s",GetName());
-       if (fTrH1 && fHits2) {
-           branch = fTrH1->GetBranch(branchname);
-           if (branch) branch->SetAddress(&fHits2);
-       }
-       if (fTrH1 && fPadHits2) {
-           branch = fTrH1->GetBranch("MUONCluster");
-           if (branch) branch->SetAddress(&fPadHits2);
-       }
-// test
-       //Int_t ntracks1 =(Int_t)fTrH1->GetEntries();
-       //printf("background - ntracks1 - %d\n",ntracks1);
-    }
-    //
-    // loop over cathodes
-    //
-    AliMUONHitMap* hm;
-    Int_t countadr=0;
-    for (int icat=0; icat<2; icat++) { 
-       Int_t counter=0;
-       Int_t * nmuon = new Int_t [AliMUONConstants::NCh()];
-       for (Int_t i =0; i<AliMUONConstants::NCh(); i++) {
-           iChamber=(AliMUONChamber*) (*fChambers)[i];
-           if (iChamber->Nsec()==1 && icat==1) {
-               continue;
-           } else {
-               segmentation=iChamber->SegmentationModel(icat+1);
-           }
-           hitMap[i] = new AliMUONHitMapA1(segmentation, list);
-           nmuon[i]=0;
-       }
-       //printf("Start loop over tracks \n");     
-//
-//   Loop over tracks
-//
-
-       TTree *treeH = gAlice->TreeH();
-       Int_t ntracks =(Int_t) treeH->GetEntries();
-       Int_t jj;
-       Float_t ** xhit = new Float_t * [2];
-       for (jj=0; jj<2; jj++) xhit[jj] = new Float_t  [AliMUONConstants::NCh()];
-       Float_t ** yhit = new Float_t * [2];
-       for (jj=0; jj<2; jj++) yhit[jj] = new Float_t  [AliMUONConstants::NCh()];
-       
-       for (Int_t track=0; track<ntracks; track++) {
-           gAlice->ResetHits();
-           treeH->GetEvent(track);
-           
-//
-//   Loop over hits
-           for(AliMUONHit* mHit=(AliMUONHit*)pMUON->FirstHit(-1); 
-               mHit;
-               mHit=(AliMUONHit*)pMUON->NextHit()) 
-           {
-               Int_t   nch   = mHit->fChamber-1;  // chamber number
-               if (nch > AliMUONConstants::NCh()-1) continue;
-               iChamber = &(pMUON->Chamber(nch));
-                // new 17.07.99
-               if (addBackground) {
-
-                   if (mHit->fParticle == kMuonPlus 
-                       || mHit->fParticle == kMuonMinus) {
-                       xhit[nch][nmuon[nch]]=mHit->fX;
-                       yhit[nch][nmuon[nch]]=mHit->fY;
-                       nmuon[nch]++;
-                       if (nmuon[nch] >2) printf("nmuon %d\n",nmuon[nch]);
-                   }
-               }
-               
-
-
-               
-//
-// Loop over pad hits
-               for (AliMUONPadHit* mPad=
-                        (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHits);
-                    mPad;
-                    mPad=(AliMUONPadHit*)pMUON->NextPad(fPadHits))
-               {
-                   Int_t cathode  = mPad->fCathode;    // cathode number
-                   Int_t ipx      = mPad->fPadX;       // pad number on X
-                   Int_t ipy      = mPad->fPadY;       // pad number on Y
-                   Int_t iqpad    = Int_t(mPad->fQpad);// charge per pad
-//
-//
-                   
-                   if (cathode != (icat+1)) continue;
-                   // fill the info array
-                   Float_t thex, they;
-                   segmentation=iChamber->SegmentationModel(cathode);
-                   segmentation->GetPadCxy(ipx,ipy,thex,they);
-//                 Float_t rpad=TMath::Sqrt(thex*thex+they*they);
-//                 if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
-
-                   new((*pAddress)[countadr++]) TVector(2);
-                   TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
-                   trinfo(0)=(Float_t)track;
-                   trinfo(1)=(Float_t)iqpad;
-
-                   digits[0]=ipx;
-                   digits[1]=ipy;
-                   digits[2]=iqpad;
-                   digits[3]=iqpad;
-                   if (mHit->fParticle == kMuonPlus ||
-                       mHit->fParticle == kMuonMinus) {
-                       digits[4]=mPad->fHitNumber;
-                   } else digits[4]=-1;
-
-                   AliMUONTransientDigit* pdigit;
-                   // build the list of fired pads and update the info
-                   if (!hitMap[nch]->TestHit(ipx, ipy)) {
-
-                       list->AddAtAndExpand(
-                           new AliMUONTransientDigit(nch,digits),counter);
-                       
-                       hitMap[nch]->SetHit(ipx, ipy, counter);
-                       counter++;
-                       pdigit=(AliMUONTransientDigit*)list->At(list->GetLast());
-                       // list of tracks
-                       TObjArray *trlist=(TObjArray*)pdigit->TrackList();
-                       trlist->Add(&trinfo);
-                   } else {
-                       pdigit=(AliMUONTransientDigit*) hitMap[nch]->GetHit(ipx, ipy);
-                       // update charge
-                       (*pdigit).fSignal+=iqpad;
-                       (*pdigit).fPhysics+=iqpad;                      
-                       // update list of tracks
-                       TObjArray* trlist=(TObjArray*)pdigit->TrackList();
-                       Int_t lastEntry=trlist->GetLast();
-                       TVector *pTrack=(TVector*)trlist->At(lastEntry);
-                       TVector &ptrk=*pTrack;
-                       Int_t lastTrack=Int_t(ptrk(0));
-                       Int_t lastCharge=Int_t(ptrk(1));
-                       if (lastTrack==track) {
-                           lastCharge+=iqpad;
-                           trlist->RemoveAt(lastEntry);
-                           trinfo(0)=lastTrack;
-                           trinfo(1)=lastCharge;
-                           trlist->AddAt(&trinfo,lastEntry);
-                       } else {
-                           trlist->Add(&trinfo);
-                       }
-                       // check the track list
-                       Int_t nptracks=trlist->GetEntriesFast();
-                       if (nptracks > 2) {
-                           for (Int_t tr=0;tr<nptracks;tr++) {
-                               TVector *ppTrack=(TVector*)trlist->At(tr);
-                               TVector &pptrk=*ppTrack;
-                               trk[tr]=Int_t(pptrk(0));
-                               chtrk[tr]=Int_t(pptrk(1));
-                           }
-                       } // end if nptracks
-                   } //  end if pdigit
-               } //end loop over clusters
-           } // hit loop
-       } // track loop
-
-       // open the file with background
-       
-       if (addBackground) {
-           ntracks =(Int_t)fTrH1->GetEntries();
-//
-//   Loop over tracks
-//
-           for (Int_t track=0; track<ntracks; track++) {
-
-               if (fHits2)       fHits2->Clear();
-               if (fPadHits2)   fPadHits2->Clear();
-
-               fTrH1->GetEvent(track);
-//
-//   Loop over hits
-               AliMUONHit* mHit;
-               for(int i=0;i<fHits2->GetEntriesFast();++i) 
-               {       
-                   mHit=(AliMUONHit*) (*fHits2)[i];
-                   Int_t   nch   = mHit->fChamber-1;  // chamber number
-                   if (nch >9) continue;
-                   iChamber = &(pMUON->Chamber(nch));
-                   Int_t rmin = (Int_t)iChamber->RInner();
-                   Int_t rmax = (Int_t)iChamber->ROuter();
-                    Float_t xbgr=mHit->fX;
-                   Float_t ybgr=mHit->fY;
-                   Bool_t cond=kFALSE;
-                   
-                   for (Int_t imuon =0; imuon < nmuon[nch]; imuon++) {
-                       Float_t dist= (xbgr-xhit[nch][imuon])*(xbgr-xhit[nch][imuon])
-                           +(ybgr-yhit[nch][imuon])*(ybgr-yhit[nch][imuon]);
-                       if (dist<100) cond=kTRUE;
-                   }
-                   if (!cond) continue;
-                   
-//
-// Loop over pad hits
-                   for (AliMUONPadHit* mPad=
-                            (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHits2);
-                        mPad;
-                        mPad=(AliMUONPadHit*)pMUON->NextPad(fPadHits2))
-                   {
-                       //                  mPad = (AliMUONPadHit*) (*fPadHits2)[j];
-                       Int_t cathode  = mPad->fCathode;    // cathode number
-                       Int_t ipx      = mPad->fPadX;       // pad number on X
-                       Int_t ipy      = mPad->fPadY;       // pad number on Y
-                       Int_t iqpad    = Int_t(mPad->fQpad);// charge per pad
-
-                       if (cathode != (icat+1)) continue;
-                       Float_t thex, they;
-                       segmentation=iChamber->SegmentationModel(cathode);
-                       segmentation->GetPadCxy(ipx,ipy,thex,they);
-                       Float_t rpad=TMath::Sqrt(thex*thex+they*they);
-                       if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
-                       new((*pAddress)[countadr++]) TVector(2);
-                       TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
-                       trinfo(0)=-1;  // tag background
-                       trinfo(1)=-1;
-                       
-                       digits[0]=ipx;
-                       digits[1]=ipy;
-                       digits[2]=iqpad;
-                       digits[3]=0;
-                       digits[4]=-1;
-                       
-                       AliMUONTransientDigit* pdigit;
-                       // build the list of fired pads and update the info
-                       if (!hitMap[nch]->TestHit(ipx, ipy)) {
-                           list->AddAtAndExpand(new AliMUONTransientDigit(nch,digits),counter);
-                           
-                           hitMap[nch]->SetHit(ipx, ipy, counter);
-                           counter++;
-                           
-                           pdigit=(AliMUONTransientDigit*)list->At(list->GetLast());
-                           // list of tracks
-                           TObjArray *trlist=(TObjArray*)pdigit->
-                               TrackList();
-                           trlist->Add(&trinfo);
-                       } else {
-                           pdigit=(AliMUONTransientDigit*) hitMap[nch]->GetHit(ipx, ipy);
-                           // update charge
-                           (*pdigit).fSignal+=iqpad;
-                           
-                           // update list of tracks
-                           TObjArray* trlist=(TObjArray*)pdigit->
-                               TrackList();
-                           Int_t lastEntry=trlist->GetLast();
-                           TVector *pTrack=(TVector*)trlist->
-                               At(lastEntry);
-                           TVector &ptrk=*pTrack;
-                           Int_t lastTrack=Int_t(ptrk(0));
-                           if (lastTrack==-1) {
-                               continue;
-                           } else {
-                               trlist->Add(&trinfo);
-                           }
-                               // check the track list
-                           Int_t nptracks=trlist->GetEntriesFast();
-                           if (nptracks > 0) {
-                               for (Int_t tr=0;tr<nptracks;tr++) {
-                                   TVector *ppTrack=(TVector*)trlist->At(tr);
-                                   TVector &pptrk=*ppTrack;
-                                   trk[tr]=Int_t(pptrk(0));
-                                   chtrk[tr]=Int_t(pptrk(1));
-                               }
-                           } // end if nptracks
-                       } //  end if pdigit
-                   } //end loop over clusters
-               } // hit loop
-           } // track loop
-           //Int_t nentr2=list->GetEntriesFast();
-           //printf(" \n counter2, nentr2 %d %d \n",counter,nentr2);
-           TTree *fAli=gAlice->TreeK();
-            TFile *file=NULL;
-           
-           if (fAli) file =fAli->GetCurrentFile();
-           file->cd();
-       } // if addBackground
-       delete [] xhit;
-       delete [] yhit;
-
-       Int_t tracks[10];
-       Int_t charges[10];
-       Int_t nentries=list->GetEntriesFast();
-       
-       for (Int_t nent=0;nent<nentries;nent++) {
-           AliMUONTransientDigit *address=(AliMUONTransientDigit*)list->At(nent);
-           if (address==0) continue; 
-           Int_t ich=address->fChamber;
-           Int_t q=address->fSignal; 
-           iChamber=(AliMUONChamber*) (*fChambers)[ich];
-//
-//  Digit Response (noise, threshold, saturation, ...)
-//             if (address->fPhysics !=0 ) address->fPhysics+=(Int_t)Noise; 
-           AliMUONResponse * response=iChamber->ResponseModel();
-           q=response->DigitResponse(q);
-           
-           if (!q) continue;
-           
-           digits[0]=address->fPadX;
-           digits[1]=address->fPadY;
-           digits[2]=q;
-           digits[3]=address->fPhysics;
-           digits[4]=address->fHit;
-           
-           TObjArray* trlist=(TObjArray*)address->TrackList();
-           Int_t nptracks=trlist->GetEntriesFast();
-           //printf("nptracks, trlist   %d  %p\n",nptracks,trlist);
-
-           // this was changed to accomodate the real number of tracks
-           if (nptracks > 10) {
-               cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
-               nptracks=10;
-           }
-           if (nptracks > 2) {
-               printf("Attention - nptracks > 2  %d \n",nptracks);
-               printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
-           }
-           for (Int_t tr=0;tr<nptracks;tr++) {
-               TVector *ppP=(TVector*)trlist->At(tr);
-               if(!ppP ) printf("ppP - %p\n",ppP);
-               TVector &pp  =*ppP;
-               tracks[tr]=Int_t(pp(0));
-               charges[tr]=Int_t(pp(1));
-                //printf("tracks, charges - %d %d\n",tracks[tr],charges[tr]);
-           }      //end loop over list of tracks for one pad
-            // Sort list of tracks according to charge
-           if (nptracks > 1) {
-               SortTracks(tracks,charges,nptracks);
-           }
-           if (nptracks < 10 ) {
-               for (Int_t i=nptracks; i<10; i++) {
-                   tracks[i]=0;
-                   charges[i]=0;
-               }
-           }
-           
-           // fill digits
-           pMUON->AddDigits(ich,tracks,charges,digits);
-           // delete trlist;
-       }
-       //cout<<"I'm out of the loops for digitisation"<<endl;
-       //      gAlice->GetEvent(nev);
-       gAlice->TreeD()->Fill();
-       pMUON->ResetDigits();
-       list->Delete();
-
-       
-       for(Int_t ii=0;ii<AliMUONConstants::NCh();++ii) {
-           if (hitMap[ii]) {
-               hm=hitMap[ii];
-               delete hm;
-               hitMap[ii]=0;
-           }
-       }
-       delete [] nmuon;    
-    } //end loop over cathodes
-    delete [] hitMap;
-    char hname[30];
-    sprintf(hname,"TreeD%d",nev);
-    gAlice->TreeD()->Write(hname);
-    // reset tree
-    gAlice->TreeD()->Reset();
-    delete list;
-    
-    pAddress->Delete();
-    // gObjectTable->Print();
-}
-
-void AliMUON::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
-{
-  //
-  // Sort the list of tracks contributing to a given digit
-  // Only the 3 most significant tracks are acctually sorted
-  //
-  
-  //
-  //  Loop over signals, only 3 times
-  //
-  
-  Int_t qmax;
-  Int_t jmax;
-  Int_t idx[3] = {-2,-2,-2};
-  Int_t jch[3] = {-2,-2,-2};
-  Int_t jtr[3] = {-2,-2,-2};
-  Int_t i,j,imax;
-  
-  if (ntr<3) imax=ntr;
-  else imax=3;
-  for(i=0;i<imax;i++){
-    qmax=0;
-    jmax=0;
-    
-    for(j=0;j<ntr;j++){
-      
-      if((i == 1 && j == idx[i-1]) 
-        ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
-      
-      if(charges[j] > qmax) {
-       qmax = charges[j];
-       jmax=j;
-      }       
-    } 
-    
-    if(qmax > 0) {
-      idx[i]=jmax;
-      jch[i]=charges[jmax]; 
-      jtr[i]=tracks[jmax]; 
-    }
-    
-  } 
-  
-  for(i=0;i<3;i++){
-    if (jtr[i] == -2) {
-         charges[i]=0;
-         tracks[i]=0;
-    } else {
-         charges[i]=jch[i];
-         tracks[i]=jtr[i];
-    }
-  }
-
-}
-
 //___________________________________________
 void AliMUON::Trigger(Int_t nev){
 // call the Trigger Algorithm and fill TreeR
@@ -1371,15 +967,31 @@ void AliMUON::Trigger(Int_t nev){
   ResetTrigger();
   char hname[30];
   sprintf(hname,"TreeR%d",nev);
-  gAlice->TreeR()->Write(hname);
+  gAlice->TreeR()->Write(hname,TObject::kOverwrite);
   gAlice->TreeR()->Reset();
   printf("\n End of trigger for event %d", nev);
 }
 
 
 //____________________________________________
-void AliMUON::FindClusters(Int_t nev,Int_t lastEntry)
+void AliMUON::Digits2Reco()
 {
+  FindClusters();
+  Int_t nev = gAlice->GetHeader()->GetEvent();
+  gAlice->TreeR()->Fill();
+  char hname[30];
+  sprintf(hname,"TreeR%d", nev);
+  gAlice->TreeR()->Write(hname);
+  gAlice->TreeR()->Reset();
+  ResetRawClusters();        
+  printf("\n End of cluster finding for event %d", nev);
+}
+
+void AliMUON::FindClusters()
+{
+//
+//  Perform cluster finding
+//
     TClonesArray *dig1, *dig2;
     Int_t ndig, k;
     dig1 = new TClonesArray("AliMUONDigit",1000);
@@ -1388,24 +1000,25 @@ void AliMUON::FindClusters(Int_t nev,Int_t lastEntry)
 //
 // Loop on chambers and on cathode planes
 //
+    ResetRawClusters();        
+    for (Int_t ich = 0; ich < 10; ich++) {
+       AliMUONChamber* iChamber = (AliMUONChamber*) (*fChambers)[ich];
+       AliMUONClusterFinderVS* rec = iChamber->ReconstructionModel();
     
-    for (Int_t ich=0;ich<10;ich++) {
-       AliMUONChamber* iChamber=(AliMUONChamber*) (*fChambers)[ich];
-       AliMUONClusterFinderVS* rec = iChamber->ReconstructionModel();    
        gAlice->ResetDigits();
-       gAlice->TreeD()->GetEvent(lastEntry);
-       TClonesArray *muonDigits  = this->DigitsAddress(ich);
+       gAlice->TreeD()->GetEvent(0);
+       TClonesArray *muonDigits = this->DigitsAddress(ich);
        ndig=muonDigits->GetEntriesFast();
        printf("\n 1 Found %d digits in %p %d", ndig, muonDigits,ich);
        TClonesArray &lhits1 = *dig1;
-       Int_t n=0;
-       for (k=0; k<ndig; k++) {
-           digit=(AliMUONDigit*) muonDigits->UncheckedAt(k);
-           if (rec->TestTrack(digit->fTracks[0]))
+       Int_t n = 0;
+       for (k = 0; k < ndig; k++) {
+           digit = (AliMUONDigit*) muonDigits->UncheckedAt(k);
+           if (rec->TestTrack(digit->Track(0)))
                new(lhits1[n++]) AliMUONDigit(*digit);
        }
        gAlice->ResetDigits();
-       gAlice->TreeD()->GetEvent(lastEntry+1);
+       gAlice->TreeD()->GetEvent(1);
        muonDigits  = this->DigitsAddress(ich);
        ndig=muonDigits->GetEntriesFast();
        printf("\n 2 Found %d digits in %p %d", ndig, muonDigits, ich);
@@ -1414,7 +1027,7 @@ void AliMUON::FindClusters(Int_t nev,Int_t lastEntry)
        
        for (k=0; k<ndig; k++) {
            digit= (AliMUONDigit*) muonDigits->UncheckedAt(k);
-           if (rec->TestTrack(digit->fTracks[0]))
+           if (rec->TestTrack(digit->Track(0)))
            new(lhits2[n++]) AliMUONDigit(*digit);
        }
 
@@ -1425,29 +1038,20 @@ void AliMUON::FindClusters(Int_t nev,Int_t lastEntry)
        dig1->Delete();
        dig2->Delete();
     } // for ich
-    gAlice->TreeR()->Fill();
-    ResetRawClusters();
-    char hname[30];
-    sprintf(hname,"TreeR%d",nev);
-    gAlice->TreeR()->Write(hname);
-    gAlice->TreeR()->Reset();
-    printf("\n End of cluster finding for event %d", nev);
-    
     delete dig1;
     delete dig2;
-    //gObjectTable->Print();
 }
  
-
+#ifdef never
 void AliMUON::Streamer(TBuffer &R__b)
 {
    // Stream an object of class AliMUON.
-      AliMUONChamber       *iChamber;
+      AliMUONChamber        *iChamber;
       AliMUONTriggerCircuit *iTriggerCircuit;
-      AliMUONSegmentation  *segmentation;
-      AliMUONResponse      *response;
-      TClonesArray         *digitsaddress;
-      TClonesArray         *rawcladdress;
+      AliSegmentation       *segmentation;
+      AliMUONResponse       *response;
+      TClonesArray          *digitsaddress;
+      TClonesArray          *rawcladdress;
       Int_t i;
       if (R__b.IsReading()) {
          Version_t R__v = R__b.ReadVersion(); if (R__v) { }
@@ -1549,6 +1153,8 @@ void AliMUON::Streamer(TBuffer &R__b)
          }
       }
 }
+#endif
+
 AliMUONPadHit* AliMUON::FirstPad(AliMUONHit*  hit, TClonesArray *clusters) 
 {
 //
@@ -1556,9 +1162,9 @@ AliMUONPadHit* AliMUON::FirstPad(AliMUONHit*  hit, TClonesArray *clusters)
     // Return the address of the first padhit for hit
     TClonesArray *theClusters = clusters;
     Int_t nclust = theClusters->GetEntriesFast();
-    if (nclust && hit->fPHlast > 0) {
-       AliMUON::fMaxIterPad=hit->fPHlast;
-       AliMUON::fCurIterPad=hit->fPHfirst;
+    if (nclust && hit->PHlast() > 0) {
+       AliMUON::fMaxIterPad=hit->PHlast();
+       AliMUON::fCurIterPad=hit->PHfirst();
        return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
     } else {
        return 0;
@@ -1567,6 +1173,8 @@ AliMUONPadHit* AliMUON::FirstPad(AliMUONHit*  hit, TClonesArray *clusters)
 
 AliMUONPadHit* AliMUON::NextPad(TClonesArray *clusters) 
 {
+// Get next pad (in iterator) 
+//
     AliMUON::fCurIterPad++;
     if (AliMUON::fCurIterPad <= AliMUON::fMaxIterPad) {
        return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
@@ -1578,6 +1186,9 @@ AliMUONPadHit* AliMUON::NextPad(TClonesArray *clusters)
 
 AliMUONRawCluster *AliMUON::RawCluster(Int_t ichamber, Int_t icathod, Int_t icluster)
 {
+//
+//  Return rawcluster (icluster) for chamber ichamber and cathode icathod
+//  Obsolete ??
     TClonesArray *muonRawCluster  = RawClustAddress(ichamber);
     ResetRawClusters();
     TTree *treeR = gAlice->TreeR();
@@ -1591,6 +1202,20 @@ AliMUONRawCluster *AliMUON::RawCluster(Int_t ichamber, Int_t icathod, Int_t iclu
     
     return  mRaw;
 }
+void   AliMUON::SetMerger(AliMUONMerger* merger)
+{
+// Set pointer to merger 
+    fMerger = merger;
+}
+
+AliMUONMerger*  AliMUON::Merger()
+{
+// Return pointer to merger
+    return fMerger;
+}
+
+
 
 AliMUON& AliMUON::operator = (const AliMUON& rhs)
 {