void MUONcombi (Int_t evNumber=0) 
{
/////////////////////////////////////////////////////////////////////////
//   This macro is a small example of a ROOT macro
//   illustrating how to read the output of GALICE
//   and fill some histograms.
//   
//     Root > .L anal.C   //this loads the macro in memory
//     Root > anal();     //by default process first event   
//     Root > anal(2);    //process third event
//Begin_Html
/*
<img src="gif/anal.gif">
*/
//End_Html
/////////////////////////////////////////////////////////////////////////
// Dynamically link some shared libs
   if (gClassTable->GetID("AliRun") < 0) {
      gROOT->LoadMacro("loadlibs.C");
      loadlibs();
   } else {
      delete gAlice;
      gAlice = 0;
   }
//
// Connect the Root Galice file containing Geometry, Kine and Hits
   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
   if (!file) file = new TFile("galice.root");

// Get AliRun object from file or create it if not on file
   if (!gAlice) {
      gAlice = (AliRun*)file->Get("gAlice");
      if (gAlice) printf("AliRun object found on file\n");
      if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
   }
   TClonesArray *PartArray = gAlice->Particles();
   TParticle *Part;
   

// Import the Kine and Hits Trees for the event evNumber in the file
   Int_t nparticles = gAlice->GetEvent(evNumber);
   if (nparticles <= 0) return;
//
//
   TH1F *dmass = new TH1F("dmass","Dimuon-Mass Distribution"
			 ,25,0.,5.);
   TH1F *dmassc = new TH1F("dmassc","Dimuon-Mass Distribution"
			 ,50,0.,10.);
   TH1F *dmassd = new TH1F("dmassd","Dimuon-Mass Distribution"
			 ,50,0.,10.);

   TH1F *pt    = new TH1F("pt","pT-single"
			 ,50,0.,10.);
//
//  Generator Loop
//
   
   AliGenCocktailEntry *Entry, *e1, *e2;
   AliGenCocktail*  Cocktail = (AliGenCocktail*) gAlice->Generator();
//                   Single Generator
   for (Entry=Cocktail->FirstGenerator();
	Entry;
	Entry=Cocktail->NextGenerator()
       ) {
       Entry->PrintInfo();
     }
//                   Pairs of Generators
//
// Initialize Combinator 
   AliDimuCombinator* Combinator = new AliDimuCombinator(PartArray);
   Combinator->SetEtaCut(-5, 5);
   Combinator->SetPtMin(0.0);   

   Int_t i=0;
//
// Single Muon  Loop  
//
   Combinator->ResetRange();
   


   for(Muon=Combinator->FirstMuon();
       Muon; 
       Muon=Combinator->NextMuon()) {
//
       Int_t chfirst= Muon->GetFirstDaughter();
       Int_t chlast = Muon->GetLastDaughter();
       Int_t parent = Muon->GetFirstMother();
       Float_t ptm  = Muon->Pt();
       Float_t eta  = Muon->Eta();
       
       printf("\n Particle %d Parent %d first child %d last child %d",
	      i,parent, chfirst, chlast);
       printf("\n Particle pt, eta: %f , %f ",pt,eta);
       i++;       
   }
//
// Di-Muon Loop

   Float_t pt1,pt2;
   TParticle* Muon1, *Muon2;

   Combinator->ResetRange();
/*
   for (Combinator->FirstMuonPairSelected(Muon1,Muon2);
	(Muon1 && Muon2);
	Combinator->NextMuonPairSelected(Muon1,Muon2))
   {
       pt1=Muon1->Pt();
       pt2=Muon2->Pt();       
       Float_t mass=Combinator->Mass(Muon1, Muon2);
       Float_t wgt =Combinator->Weight(Muon1, Muon2);
       pt->Fill(pt1, wgt);
       pt->Fill(pt2, wgt);
       Float_t smeared_mass=mass;
       Combinator->SmearGauss(0.05*mass, smeared_mass);
       if (Combinator->Correlated(Muon1, Muon2)) {
	   dmassc->Fill(mass, wgt);
       } else {
	   dmass->Fill(mass, wgt);
       }
    }
*/
//
//   Dimuon Loop controlled by Generator Loop
//
   for (Cocktail->FirstGeneratorPair(e1,e2);
	(e1&&e2);
	Cocktail->NextGeneratorPair(e1,e2)
       ){
       printf("\n ----------------------------------------------------");
       e1->PrintInfo();
       e2->PrintInfo();
       printf("\n ----------------------------------------------------");
       Combinator->SetFirstRange (e1->GetFirst(), e1->GetLast());
       Combinator->SetSecondRange(e2->GetFirst(), e2->GetLast());
       Combinator->SetRate(e1->Rate(), e2->Rate());
       for (Combinator->FirstMuonPairSelected(Muon1,Muon2);
	    (Muon1 && Muon2);
	    Combinator->NextMuonPairSelected(Muon1,Muon2))
       {
	   pt1=Muon1->Pt();
	   pt2=Muon2->Pt();       
	   Float_t mass=Combinator->Mass(Muon1, Muon2);
	   Float_t wgt =Combinator->Weight(Muon1, Muon2);
	   pt->Fill(pt1, wgt);
	   pt->Fill(pt2, wgt);
	   Float_t smeared_mass=mass;
	   Combinator->SmearGauss(0.05*mass, smeared_mass);
	   Float_t DecayWeight=
	       Combinator->Decay_Prob(Muon2)*
	       Combinator->Decay_Prob(Muon2);
//	   if (TMath::Min(pt1,pt2) > -0.5*TMath::Max(pt1,pt2)+2.) {
	       if (Combinator->Correlated(Muon1, Muon2)) {
		   dmassc->Fill(smeared_mass, wgt);
	       }
		   dmass->Fill(smeared_mass, wgt);
		   dmassd->Fill(smeared_mass, wgt*DecayWeight);
//	   }
       } // Dimuon Loop
   }// Generator Loop
//
//Create a canvas, set the view range, show histograms
//
   TCanvas *c1 = new TCanvas("c1","Dimuon Plots",400,10,600,700);
   TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
   TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
   TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
   TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
   pad11->SetFillColor(11);
   pad12->SetFillColor(11);
   pad13->SetFillColor(11);
   pad14->SetFillColor(11);
   pad11->Draw();
   pad12->Draw();
   pad13->Draw();
   pad14->Draw();
   
   pad11->cd();
   pt->SetFillColor(42);
   pt->SetXTitle("pT (GeV)");
   pt->Draw();

   pad12->cd();
   pad12->SetLogy(0);
   dmass->SetFillColor(42);
   dmass->SetXTitle("m (GeV)");
   dmass->Draw();

   pad13->cd();
   pad13->SetLogy(0);
   dmassc->SetFillColor(42);
   dmassc->SetXTitle("m (GeV)");
   dmassc->Draw();

   pad14->cd();
   pad14->SetLogy(0);
   dmassd->SetFillColor(42);
   dmassd->SetXTitle("m (GeV)");
   dmassd->Draw();


}