void toy(){


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

   TLorentzVector W(0.0, 0.0, 0.0, 4.42);

   //(Momentum, Energy units are Gev/C, GeV)
   Double_t masses[3] = { 3.686, 0.139, 0.139} ;

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

  double x_min = (masses[0]+masses[1])*(masses[0]+masses[1]);
  double x_max = (4.42-masses[1])*(4.42-masses[1]);

  TFile* output = new TFile("toymc.root", "recreate");

   TH2F *h2 = new TH2F("h2","h2", 100,x_min,x_max, 100,x_min,x_max);

   for (Int_t n=0;n<100000;n++) {
      Double_t weight = event.Generate();

      TLorentzVector *Ppsip= event.GetDecay(0);

      TLorentzVector *pPip    = event.GetDecay(1);
      TLorentzVector *pPim    = event.GetDecay(2);

      TLorentzVector pPpsipPip = *Ppsip + *pPip;
      TLorentzVector pPpsipPim = *Ppsip + *pPim;

      h2->Fill(pPpsipPip.M2() ,pPpsipPim.M2() ,weight);
   }

   output->Write();
   output->Close();
}