// c++ include
#include <iostream>
#include <iomanip>
#include <fstream>

// ROOT include
#include <TROOT.h>
#include <TSystem.h>
#include <TApplication.h>
#include <TFile.h>
#include <TH1.h>
#include <TNtuple.h>
#include <TGraphErrors.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TStyle.h>
#include <TObject.h>
#include <TMinuit.h>
#include <TLegend.h>
#include <TPaveText.h>
#include <TLorentzVector.h>
#include <TGenPhaseSpace.h>


void PhaseSpace(int events = 10000000)
{
    
  if (!gROOT->GetClass("TGenPhaseSpace")) gSystem->Load("libPhysics");
  
  TCanvas *c2 = new TCanvas("c2","c2",1,1,400,400);
  c2->SetFillColor(10);
  
   
  // PHYSICS
  Double_t mmu = 0.105658367;
  
  TLorentzVector *Pmup = new TLorentzVector(0.0, 0.0, 0.0, 10*mmu);

  //(Momentum, Energy units are Gev/C, GeV)
  Double_t masses[3] = {0.0, 0.0, 0.0} ;
  
  // Generate Phass space for e+,e-,e+,nu_e,nub_mu
  TGenPhaseSpace event;
  event.SetDecay((*Pmup), 3, masses);
  
  TH1F *V      = new TH1F("V", "V", 1,0.5,1.5);

  // BLHA interface
  double* p_momenta = new double[20];

 
  for (Int_t n=0;n<events;n++) {
    
    Double_t weight = event.Generate();
    
    TLorentzVector *Pep1   = event.GetDecay(0);
    TLorentzVector *Pnue   = event.GetDecay(1);
    TLorentzVector *Pnubmu = event.GetDecay(2);
    
    V -> Fill(1.0, weight); 
    
  }

  c2->cd();
  V->Scale(1.0/V->GetEntries());
  std::cout<<"Volume="<<V->Integral()<<std::endl;
  V->Draw();

}