void f22806(UInt_t NrOfEvents = 1000000) {
   if (!gROOT->GetClass("TGenPhaseSpace")) gSystem->Load("libPhysics");

   TLorentzVector target(0.0, 0.0, 0.0, 1.879);
   TLorentzVector beam(0.0, 0.0, 2.2, 2.2);
   TLorentzVector W = beam + target;

   //(Momentum, Energy units are Gev/C, GeV)
   Double_t masses[4] = {0.938, 0.493, 0.493, 0.938} ;

   //int detectn = 100000; //detect n will eventually be an output based on detector efficiency and angle p
   //void doit(UInt_t NrOfEvents = 1000000) //***Number of events can be changed here or when you want to run by "doit(...)"
   //{
   int c=0;
   TRandom3 p(-1);
   for (UInt_t iEvent = 0; iEvent < NrOfEvents; ++iEvent) {
      Double_t num = p.Uniform(0,360);
      if (num < 0.05){
         c=c+1;
      }
   }
   cout << "c="<< c << endl;
}