#include "TLorentzVector.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TGenPhaseSpace.h"
#include "TRandom3.h"

#include <iostream>

using namespace std;

void test(){
  double E=11.;
  double Mp=0.9383;
  double Mpi=0.134976;
  double Me=0.511E-3;


  TLorentzVector Pbeam(0,0,sqrt(E*E-Me*Me),E);
  TLorentzVector Ptarget(0,0,0,Mp);

  TLorentzVector P0=(Pbeam+Ptarget);

  TLorentzVector Pe,Pp,Ppi;
  TGenPhaseSpace generator;

  TH1D *ht=new TH1D("ht","ht",250,0,5);

  double dmasses[3]={Me,Mp,Mpi};

  if(!generator.SetDecay(P0,3,dmasses)){
    cout<<"The process is kinematically forbidden!"<<endl;
    exit(1);
  }

  double mW=generator.GetWtMax();
  double w,W;
  double t,Q2;
  bool flag_printError;
  gRandom->SetSeed(0);
  for (int ii=0;ii<20000000;ii++){
    w=generator.Generate();

    if ((w>mW)&&(ii<1000)){
      cout<<"Error! "<<w<<" "<<mW<<endl;
    }
    if (w<gRandom->Uniform(0,mW)){
      ii--;
      continue;
    }

    Pe=(*generator.GetDecay(0));
    Pp=(*generator.GetDecay(1));
    Ppi=(*generator.GetDecay(2));
    
    t=-(Ptarget-Pp).M2(); //magnitude of momentum transferred^2 between target and recoil
    Q2=-(Pe-Pbeam).M2();

    W=(Pp+Ppi).M2();

    
    if ((fabs(W-4)<0.5)&&(fabs(Q2-1)<0.2)){
      ht->Fill(t);
    }
  }

  ht->Draw();



}