]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUON.cxx
Possibility to investigate a primary of not yet loaded particle (I.Hrivnacova)
[u/mrichter/AliRoot.git] / MUON / AliMUON.cxx
index 1ca5dd1101a16d8ee93aba1d5fda8629dc976809..b46632ee97edced58540b2600300eda65ef6d1b8 100644 (file)
  **************************************************************************/
 /*
 $Log$
+Revision 1.54  2001/08/30 09:52:12  hristov
+The operator[] is replaced by At() or AddAt() in case of TObjArray.
+
+Revision 1.53  2001/07/20 10:03:13  morsch
+Changes needed to work with Root 3.01 (substitute lhs [] operator). (Jiri Chudoba)
+
+Revision 1.52  2001/06/14 13:49:22  hristov
+Write a TreeD in SDigits2Digits method (needed to be compatible with alirun script)
+
+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
 
@@ -98,7 +126,7 @@ 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 AliMUONClusterFinderVS.
+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).
@@ -222,8 +250,10 @@ Log message added
 #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" 
@@ -267,6 +297,11 @@ AliMUON::AliMUON()
     fAccMin          = 0.;
     fAccMax          = 0.;   
     fAccCut          = kFALSE;
+    fMerger          = 0;
+    fFileName        = 0;
+    fTrH1            = 0;
+    fHits2           = 0;
+    fPadHits2        = 0;
 }
  
 //___________________________________________
@@ -294,7 +329,7 @@ AliMUON::AliMUON(const char *name, const char *title)
    Int_t i;
    
    for (i=0; i<AliMUONConstants::NCh() ;i++) {
-       (*fDchambers)[i] = new TClonesArray("AliMUONDigit",10000); 
+       fDchambers->AddAt(new TClonesArray("AliMUONDigit",10000),i); 
        fNdch[i]=0;
    }
 
@@ -303,7 +338,7 @@ AliMUON::AliMUON(const char *name, const char *title)
    fRawClusters = new TObjArray(AliMUONConstants::NTrackingCh());
 
    for (i=0; i<AliMUONConstants::NTrackingCh();i++) {
-       (*fRawClusters)[i] = new TClonesArray("AliMUONRawCluster",10000); 
+       fRawClusters->AddAt(new TClonesArray("AliMUONRawCluster",10000),i); 
        fNrawch[i]=0;
    }
 
@@ -334,12 +369,13 @@ AliMUON::AliMUON(const char *name, const char *title)
            ch = 2 * st + stCH;
 //
            if (ch < AliMUONConstants::NTrackingCh()) {
-               (*fChambers)[ch] = new AliMUONChamber(ch);
+             fChambers->AddAt(new AliMUONChamber(ch),ch);
            } else {
-               (*fChambers)[ch] = new AliMUONChamberTrigger(ch);
+             fChambers->AddAt(new AliMUONChamberTrigger(ch),ch);
            }
            
-           AliMUONChamber* chamber = (AliMUONChamber*) (*fChambers)[ch];
+        //PH       AliMUONChamber* chamber = (AliMUONChamber*) (*fChambers)[ch];
+           AliMUONChamber* chamber = (AliMUONChamber*) fChambers->At(ch);
            
            chamber->SetGid(0);
            // Default values for Z of chambers
@@ -352,6 +388,7 @@ AliMUON::AliMUON(const char *name, const char *title)
 //
        } // Chamber stCH (0, 1) in 
     }     // Station st (0...)
+//    fChambers->SetLast(AliMUONConstants::NCh());
     fMaxStepGas=0.01; 
     fMaxStepAlu=0.1; 
     fMaxDestepGas=-1;
@@ -367,10 +404,10 @@ AliMUON::AliMUON(const char *name, const char *title)
    // 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
+     fTriggerCircuits->AddAt(new AliMUONTriggerCircuit(),circ);     
 
+   }
+     fMerger = 0;
 }
  
 //___________________________________________
@@ -384,7 +421,7 @@ AliMUON::AliMUON(const AliMUON& rMUON)
 AliMUON::~AliMUON()
 {
 // Destructor
-    printf("Calling AliMUON destructor !!!\n");
+    if(fDebug) printf("%s: Calling AliMUON destructor !!!\n",ClassName());
     
     Int_t i;
     fIshunt  = 0;
@@ -454,6 +491,8 @@ AliMUON::~AliMUON()
       fTrH1->Delete();
       delete fTrH1;
     }
+
+    if (fMerger) delete fMerger;
 }
  
 //___________________________________________
@@ -475,7 +514,8 @@ void AliMUON::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
     // Add a MUON digit to the list
     //
 
-    TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
+  //PH    TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
+    TClonesArray &ldigits = *((TClonesArray*)fDchambers->At(id));
     new(ldigits[fNdch[id]++]) AliMUONDigit(tracks,charges,digits);
 }
 
@@ -486,7 +526,8 @@ void AliMUON::AddRawCluster(Int_t id, const AliMUONRawCluster& c)
     // Add a MUON digit to the list
     //
 
-    TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
+  //PH    TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
+    TClonesArray &lrawcl = *((TClonesArray*)fRawClusters->At(id));
     new(lrawcl[fNrawch[id]++]) AliMUONRawCluster(c);
 }
 
@@ -528,7 +569,7 @@ Int_t AliMUON::DistancetoPrimitive(Int_t , Int_t )
 }
 
 //___________________________________________
-void AliMUON::MakeBranch(Option_t* option, char *file)
+void AliMUON::MakeBranch(Option_t* option, const char *file)
 {
     //
     // Create Tree branches for the MUON.
@@ -539,14 +580,13 @@ void AliMUON::MakeBranch(Option_t* option, char *file)
     
     AliDetector::MakeBranch(option,file);
     
-    char *cD = strstr(option,"D");
-    char *cR = strstr(option,"R");
-    char *cH = strstr(option,"H");
+    const char *cD = strstr(option,"D");
+    const char *cR = strstr(option,"R");
+    const char *cH = strstr(option,"H");
 
     if (fPadHits   && gAlice->TreeH() && cH) {
-      gAlice->MakeBranchInTree(gAlice->TreeH(), 
-                               branchname, &fPadHits, kBufferSize, file) ;       
-         printf("Making Branch %s for clusters\n",branchname);
+      MakeBranchInTree(gAlice->TreeH(), 
+                       branchname, &fPadHits, kBufferSize, file);
     }
     
     if (cD) {
@@ -558,8 +598,8 @@ void AliMUON::MakeBranch(Option_t* option, char *file)
       for (i=0; i<AliMUONConstants::NCh() ;i++) {
            sprintf(branchname,"%sDigits%d",GetName(),i+1);     
            if (fDchambers   && gAlice->TreeD()) {
-          gAlice->MakeBranchInTree(gAlice->TreeD(), 
-                                   branchname, &((*fDchambers)[i]), kBufferSize, file) ;         
+            MakeBranchInTree(gAlice->TreeD(), 
+                             branchname, &((*fDchambers)[i]), kBufferSize, file);
              printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
         }
          }     
@@ -576,8 +616,8 @@ void AliMUON::MakeBranch(Option_t* option, char *file)
       for (i=0; i<AliMUONConstants::NTrackingCh() ;i++) {
            sprintf(branchname,"%sRawClusters%d",GetName(),i+1);        
            if (fRawClusters   && gAlice->TreeR()) {
-          gAlice->MakeBranchInTree(gAlice->TreeR(), 
-                                   branchname, &((*fRawClusters)[i]), kBufferSize, file) ;       
+              MakeBranchInTree(gAlice->TreeR(), 
+                               branchname, &((*fRawClusters)[i]), kBufferSize, file);
              printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
            }   
       }
@@ -586,8 +626,8 @@ void AliMUON::MakeBranch(Option_t* option, char *file)
       //
       sprintf(branchname,"%sGlobalTrigger",GetName());
       if (fGlobalTrigger && gAlice->TreeR()) {  
-        gAlice->MakeBranchInTree(gAlice->TreeR(), 
-                                 branchname, &fGlobalTrigger, kBufferSize, file) ;       
+          MakeBranchInTree(gAlice->TreeR(), 
+                           branchname, &fGlobalTrigger, kBufferSize, file);
            printf("Making Branch %s for Global Trigger\n",branchname);
       }
       //
@@ -595,8 +635,8 @@ void AliMUON::MakeBranch(Option_t* option, char *file)
       //  
       sprintf(branchname,"%sLocalTrigger",GetName());
       if (fLocalTrigger && gAlice->TreeR()) {  
-        gAlice->MakeBranchInTree(gAlice->TreeR(), 
-                                 branchname, &fLocalTrigger, kBufferSize, file) ;        
+          MakeBranchInTree(gAlice->TreeR(), 
+                           branchname, &fLocalTrigger, kBufferSize, file);
            printf("Making Branch %s for Local Trigger\n",branchname);
       }
    }
@@ -668,7 +708,8 @@ void AliMUON::ResetDigits()
     // Reset number of digits and the digits array for this detector
     //
     for ( int i=0;i<AliMUONConstants::NCh();i++ ) {
-       if ((*fDchambers)[i])    ((TClonesArray*)(*fDchambers)[i])->Clear();
+      //PH     if ((*fDchambers)[i])    ((TClonesArray*)(*fDchambers)[i])->Clear();
+       if ((*fDchambers)[i])    ((TClonesArray*)fDchambers->At(i))->Clear();
        if (fNdch)  fNdch[i]=0;
     }
 }
@@ -679,7 +720,8 @@ void AliMUON::ResetRawClusters()
     // Reset number of raw clusters and the raw clust array for this detector
     //
     for ( int i=0;i<AliMUONConstants::NTrackingCh();i++ ) {
-       if ((*fRawClusters)[i])    ((TClonesArray*)(*fRawClusters)[i])->Clear();
+      //PH     if ((*fRawClusters)[i])    ((TClonesArray*)(*fRawClusters)[i])->Clear();
+       if ((*fRawClusters)[i])    ((TClonesArray*)fRawClusters->At(i))->Clear();
        if (fNrawch)  fNrawch[i]=0;
     }
 }
@@ -699,8 +741,10 @@ 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);
+    //PH    ((AliMUONChamber*) (*fChambers)[i])  ->SetPadSize(isec,p1,p2);
+    //PH    ((AliMUONChamber*) (*fChambers)[i+1])->SetPadSize(isec,p1,p2);
+    ((AliMUONChamber*) fChambers->At(i))  ->SetPadSize(isec,p1,p2);
+    ((AliMUONChamber*) fChambers->At(i+1))->SetPadSize(isec,p1,p2);
 }
 
 //___________________________________________
@@ -709,7 +753,8 @@ void AliMUON::SetChambersZ(const Float_t *Z)
   // Set Z values for all chambers (tracking and trigger)
   // from the array pointed to by "Z"
     for (Int_t ch = 0; ch < AliMUONConstants::NCh(); ch++)
-       ((AliMUONChamber*) ((*fChambers)[ch]))->SetZ(Z[ch]);
+      //PH     ((AliMUONChamber*) ((*fChambers)[ch]))->SetZ(Z[ch]);
+       ((AliMUONChamber*) fChambers->At(ch))->SetZ(Z[ch]);
     return;
 }
 
@@ -727,8 +772,10 @@ 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);
+    //PH    ((AliMUONChamber*) (*fChambers)[i])->SetChargeSlope(p1);
+    //PH    ((AliMUONChamber*) (*fChambers)[i+1])->SetChargeSlope(p1);
+    ((AliMUONChamber*) fChambers->At(i))->SetChargeSlope(p1);
+    ((AliMUONChamber*) fChambers->At(i+1))->SetChargeSlope(p1);
 }
 
 //___________________________________________
@@ -736,8 +783,10 @@ 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);
+    //PH    ((AliMUONChamber*) fChambers->At(i))->SetChargeSpread(p1,p2);
+    //PH    ((AliMUONChamber*) fChambers->Ati+1])->SetChargeSpread(p1,p2);
+    ((AliMUONChamber*) fChambers->At(i))->SetChargeSpread(p1,p2);
+    ((AliMUONChamber*) fChambers->At(i+1))->SetChargeSpread(p1,p2);
 }
 
 //___________________________________________
@@ -745,8 +794,10 @@ 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);
+    //PH    ((AliMUONChamber*) (*fChambers)[i])->SetSigmaIntegration(p1);
+    //PH    ((AliMUONChamber*) (*fChambers)[i+1])->SetSigmaIntegration(p1);
+    ((AliMUONChamber*) fChambers->At(i))->SetSigmaIntegration(p1);
+    ((AliMUONChamber*) fChambers->At(i+1))->SetSigmaIntegration(p1);
 }
 
 //___________________________________________
@@ -754,8 +805,10 @@ 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);
+    //PH    ((AliMUONChamber*) (*fChambers)[i])->SetMaxAdc(p1);
+    //PH    ((AliMUONChamber*) (*fChambers)[i+1])->SetMaxAdc(p1);
+    ((AliMUONChamber*) fChambers->At(i))->SetMaxAdc(p1);
+    ((AliMUONChamber*) fChambers->At(i+1))->SetMaxAdc(p1);
 }
 
 //___________________________________________
@@ -809,48 +862,51 @@ void AliMUON::SetAcceptance(Bool_t acc, Float_t angmin, Float_t angmax)
 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);
+  //PH    ((AliMUONChamber*) (*fChambers)[id])->SetSegmentationModel(isec, segmentation);
+    ((AliMUONChamber*) fChambers->At(id))->SetSegmentationModel(isec, segmentation);
 
 }
 //___________________________________________
 void   AliMUON::SetResponseModel(Int_t id, AliMUONResponse *response)
 {
 // Set the response for chamber id
-    ((AliMUONChamber*) (*fChambers)[id])->SetResponseModel(response);
+  //PH    ((AliMUONChamber*) (*fChambers)[id])->SetResponseModel(response);
+    ((AliMUONChamber*) fChambers->At(id))->SetResponseModel(response);
 }
 
 void   AliMUON::SetReconstructionModel(Int_t id, AliMUONClusterFinderVS *reconst)
 {
 // Set ClusterFinder for chamber id
-    ((AliMUONChamber*) (*fChambers)[id])->SetReconstructionModel(reconst);
+  //PH    ((AliMUONChamber*) (*fChambers)[id])->SetReconstructionModel(reconst);
+    ((AliMUONChamber*) fChambers->At(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);
+  //PH    ((AliMUONChamber*) (*fChambers)[id])->SetNsec(nsec);
+    ((AliMUONChamber*) fChambers->At(id))->SetNsec(nsec);
 }
 
 //___________________________________________
 void AliMUON::SDigits2Digits()
 {
-    Int_t evNumber2=0, ibg=0, bgr=0;
-    Int_t nbgr_ev= 0, nev=0;
-
-    Int_t nparticles = gAlice->GetNtrack();
-    cout << "nev         " <<nev<<endl;
-    cout << "nparticles  " <<nparticles<<endl; 
-    if (nparticles <= 0) return;
-
-    nbgr_ev = Int_t(nev*bgr/(evNumber2+1));
-       
-    if (ibg) {
-         printf("nbgr_ev %d\n",nbgr_ev);
-         Digitise(nev,nbgr_ev,"Add"," ","galice_bgr.root");
-    } 
-       else {
-         Digitise(nev,nbgr_ev,"rien","",""); 
+
+// 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();
 }
 
 //___________________________________________
@@ -876,7 +932,8 @@ void AliMUON::MakePadHits(Float_t xhit,Float_t yhit, Float_t zhit,
 //    if (idvol == 6) printf("\n ->Disintegration %f %f %f", xhit, yhit, eloss );
     
 
-    ((AliMUONChamber*) (*fChambers)[idvol])
+    //PH    ((AliMUONChamber*) (*fChambers)[idvol])
+    ((AliMUONChamber*) fChambers->At(idvol))
        ->DisIntegration(eloss, tof, xhit, yhit, zhit, nnew, newclust);
     Int_t ic=0;
 //    if (idvol == 6) printf("\n nnew  %d \n", nnew);
@@ -903,479 +960,6 @@ void AliMUON::MakePadHits(Float_t xhit,Float_t yhit, Float_t zhit,
     }
 }
 
-//----------------------------------------------------------------------
-
-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;
-    AliSegmentation*  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");
-    AliHitMap** hitMap= new AliHitMap* [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);
-       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);
-       }
-    }
-    //
-    // loop over cathodes
-    //
-    AliHitMap* 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;
-       Int_t nmaxmuon = 5;
-
-       Float_t ** xhit = new Float_t * [AliMUONConstants::NCh()];
-       for (jj = 0; jj < AliMUONConstants::NCh(); jj++)
-         xhit[jj] = new Float_t[nmaxmuon];
-       Float_t ** yhit = new Float_t * [AliMUONConstants::NCh()];
-       for (jj=0; jj<AliMUONConstants::NCh(); jj++)
-         yhit[jj] = new Float_t[nmaxmuon];
-
-       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->Chamber()-1;  // chamber number
-               if (nch > AliMUONConstants::NCh()-1) continue;
-//             if (nch > 9) continue;
-               iChamber = &(pMUON->Chamber(nch));
-                // new 17.07.99
-               if (addBackground) {
-
-                   if (mHit->Particle()    == kMuonPlus 
-                       || mHit->Particle() == kMuonMinus) {
-                     // mark muon hits
-                     if (nmuon[nch] < nmaxmuon) {
-                       xhit[nch][nmuon[nch]] = mHit->X();
-                       yhit[nch][nmuon[nch]] = mHit->Y();
-                       nmuon[nch]++;
-                     }
-                     // ignore muon if too many compared to nmaxmuon
-                     else printf("AliMUON::Digitize: nmuon %d ==> ignored\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->Cathode();     // cathode number
-                   Int_t ipx      = mPad->PadX();        // pad number on X
-                   Int_t ipy      = mPad->PadY();        // pad number on Y
-                   Int_t iqpad    = mPad->QPad();        // charge per pad
-//                 printf("\n Pad: %d %d %d %d", ipx, ipy, cathode,nch);
-//
-//
-                   
-                   if (cathode != (icat+1)) continue;
-                   // fill the info array
-                   segmentation = iChamber->SegmentationModel(cathode);
-                   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] = icat;
-                   digits[3] = iqpad;
-                   digits[4] = iqpad;
-
-                   if (mHit->Particle() == kMuonPlus ||
-                       mHit->Particle() == kMuonMinus) {
-                       digits[5] = mPad->HitNumber();
-                   } else digits[5] = -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).AddSignal(iqpad);
-                       (*pdigit).AddPhysicsSignal(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->Chamber()-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->X();
-                   Float_t ybgr = mHit->Y();
-                   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->Cathode();    // cathode number
-                       Int_t ipx      = mPad->PadX();       // pad number on X
-                       Int_t ipy      = mPad->PadY();       // pad number on Y
-                       Int_t iqpad    = mPad->QPad();       // charge per pad
-
-                       if (cathode != (icat+1)) 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] = icat;
-                       digits[3] = iqpad;
-                       digits[4] =  0;
-                       digits[5] = -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).AddSignal(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
-
-           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->Chamber();
-           Int_t   q = address->Signal(); 
-           iChamber = (AliMUONChamber*) (*fChambers)[ich];
-//
-//  Digit Response (noise, threshold, saturation, ...)
-
-           AliMUONResponse * response = iChamber->ResponseModel();
-           q=response->DigitResponse(q);
-           
-           if (!q) continue;
-           
-           digits[0] = address->PadX();
-           digits[1] = address->PadY();
-           digits[2] = address->Cathode();
-           digits[3] = q;
-           digits[4] = address->Physics();
-           digits[5] = address->Hit();
-           
-           TObjArray* trlist = (TObjArray*)address->TrackList();
-           Int_t nptracks = trlist->GetEntriesFast();
-
-           // 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,TObject::kOverwrite);
-    // 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
@@ -1425,7 +1009,20 @@ void AliMUON::Trigger(Int_t 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
@@ -1438,12 +1035,14 @@ 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();    
+      //PH     AliMUONChamber* iChamber = (AliMUONChamber*) (*fChambers)[ich];
+       AliMUONChamber* iChamber = (AliMUONChamber*) fChambers->At(ich);
+       AliMUONClusterFinderVS* rec = iChamber->ReconstructionModel();
+    
        gAlice->ResetDigits();
-       gAlice->TreeD()->GetEvent(lastEntry);
+       gAlice->TreeD()->GetEvent(0);
        TClonesArray *muonDigits = this->DigitsAddress(ich);
        ndig=muonDigits->GetEntriesFast();
        printf("\n 1 Found %d digits in %p %d", ndig, muonDigits,ich);
@@ -1455,7 +1054,7 @@ void AliMUON::FindClusters(Int_t nev,Int_t lastEntry)
                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);
@@ -1475,17 +1074,8 @@ 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,TObject::kOverwrite);
-    gAlice->TreeR()->Reset();
-    printf("\n End of cluster finding for event %d", nev);
-    
     delete dig1;
     delete dig2;
-    //gObjectTable->Print();
 }
  
 #ifdef never
@@ -1648,6 +1238,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)
 {
@@ -1656,21 +1260,28 @@ AliMUON& AliMUON::operator = (const AliMUON& rhs)
     return *this;
 }
 
+////////////////////////////////////////////////////////////////////////
+void AliMUON::MakeBranchInTreeD(TTree *treeD, const char *file)
+{
+    //
+    // Create TreeD branches for the MUON.
+    //
 
+  const Int_t kBufferSize = 4000;
+  char branchname[30];
+    
+  //
+  // one branch for digits per chamber
+  // 
+  for (Int_t i=0; i<AliMUONConstants::NCh() ;i++) {
+    sprintf(branchname,"%sDigits%d",GetName(),i+1);    
+    if (fDchambers && treeD) {
+      MakeBranchInTree(treeD, 
+                      branchname, &((*fDchambers)[i]), kBufferSize, file);
+      printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
+    }
+  }
+}
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+//___________________________________________