{
  gROOT.Reset();
  if (!gROOT->GetClass("TGenPhaseSpace")) gSystem.Load("libPhysics");

  // particle multiplicity
  const int Np = 5;   

  // final state particles
  double mass[Np] = { 0.938, 0.140, 0.140, 0.140, 0.140 };

  // invariant mass & 4p at hadronic CM
  double W = 1.5; 
  TLorentzVector p(0.,0.,0.,W);

  // create phase space generator and set decay
  TGenPhaseSpace gen;
  bool permitted = gen.SetDecay(p, Np, mass);
  assert(permitted);

  // max weight for specified decay
  double max_wgt = gen.GetWtMax();
  cout << "Max weight: " << max_wgt << endl;

  // fill log10(weights) histogram
  const int NDecays = 10000000;
  TH1F * hwgt = new TH1F("hwgt","phase space decay weights",500,-3,TMath::Log10(4*max_wgt));

  for(int id=0; id<NDecays; id++) {
     double w = gen.Generate();
     if(w>0) {
      hwgt->Fill(TMath::Log10(w));
     } 
 }

  // plot
  TCanvas * c = new TCanvas("c","",20,20,500,500);
  hwgt->Draw();
  hwgt->GetXaxis()->SetTitle("log10(weight)");

  TLine * limit = new TLine(TMath::Log10(max_wgt), 0, TMath::Log10(max_wgt), hwgt->GetMaximum());
  limit->SetLineColor(2);
  limit->SetLineStyle(2);
  limit->Draw();
}