//To run this script, you can use CINT or ACLIC
// -with CINT, do:
//   root > .x genie.C
//
// -with ACLIC, do
// root > gSystem->Load("libPhysics");
// root > .x genie.C+
//
#include "TGenPhaseSpace.h"
#include "TLorentzVector.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TLine.h"
#include "Riostream.h"
   
void genie() {

  // 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));

  Int_t neg=0;
  for(int id=0; id<NDecays; id++) {
     if (id%1000000 == 0) printf("id = %d\n",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();
}