 
void mc5a(int iShow) {
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);

  TH2D *hDalitzWeighted = new TH2D("hDalitzWeighted","",120,0.0,12.0,200,0.0,20);
  TH2D *hDalitzNotWeighted = new TH2D("hDalitzNotWeighted","",120,0.0,12.0,200,0.0,20);
 
  double mPro       = 0.938;
  double mPip       = 0.13957;
  double mPim       = mPip;

  TLorentzVector target(0,0,0,mPro);
  TLorentzVector beam(0.0,0.0,8.5,8.5);
  TLorentzVector pTot4Vec = beam + target;

  Double_t masses[3] = {mPro,mPip,mPim} ;

  TGenPhaseSpace event;
  event.SetDecay(pTot4Vec, 3, masses);

  int nToGen = 1000000;
  for (int i = 0; i<nToGen; i++) {
    Double_t weight = event.Generate();

    TLorentzVector *pProton = event.GetDecay(0);
    TLorentzVector *pPip    = event.GetDecay(1);
    TLorentzVector *pPim    = event.GetDecay(2);
 
    TLorentzVector pPPip   = *pProton + *pPip;
    TLorentzVector pPPim   = *pProton + *pPim;
    TLorentzVector pPipPim = *pPip + *pPim;
 
    hDalitzWeighted->Fill(pPipPim.M2() ,pPPim.M2() ,weight);
    hDalitzNotWeighted->Fill(pPipPim.M2() ,pPPim.M2());
  }
  if (iShow == 1) hDalitzWeighted->Draw("colz");
  if (iShow == 2) hDalitzWeighted->Draw("lego2");
  if (iShow == 3) hDalitzNotWeighted->Draw("lego2");
}
