diff --git a/data/emcal_pi0_0GeVto30GeV_1000kEvt.hepmc b/data/emcal_pi0_0GeVto30GeV_1000kEvt.hepmc
new file mode 100644
index 0000000000000000000000000000000000000000..112f59c27c98360178ada360580c539f2708bc3c
Binary files /dev/null and b/data/emcal_pi0_0GeVto30GeV_1000kEvt.hepmc differ
diff --git a/emcal_pi0.cxx b/emcal_pi0.cxx
index bbab569cbe9081ce92f66c013f5c1e915bc3accd..69018e7806f3cd94bd46c905b8f6b75dc3b70da0 100644
--- a/emcal_pi0.cxx
+++ b/emcal_pi0.cxx
@@ -16,7 +16,7 @@
 
 using namespace HepMC3;
 
-void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0GeVto30GeV_100kEvt.hepmc")
+void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0GeVto30GeV_1000kEvt.hepmc")
 {
   WriterAscii hepmc_output(out_fname);
   int events_parsed = 0;
@@ -25,6 +25,12 @@ void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0Ge
   // Random number generator
   TRandom *r1 = new TRandom();
 
+  // Constraining the solid angle, but larger than that subtended by the detector
+  // https://indico.bnl.gov/event/7449/contributions/35966/attachments/27177/41430/EIC-DWG-Calo-03192020.pdf
+  // See a figure on slide 26
+  double cos_theta_min = std::cos(M_PI * (120.0 / 180.0));
+  double cos_theta_max = std::cos(M_PI);
+
   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
     // FourVector(px,py,pz,e,pdgid,status)
     // type 4 is beam
@@ -37,20 +43,25 @@ void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0Ge
         FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
 
     // Define momentum
-    Double_t p = 0.0 + events_parsed * 3.e-4;
-    Double_t px;
-    Double_t py;
-    Double_t pz;
+    // Momentum starting at 0 GeV/c to 30 GeV/c
+    Double_t p        = r1->Uniform(0.0, 30.0);
+    Double_t phi      = r1->Uniform(0.0, 2.0 * M_PI);
+    Double_t costheta = r1->Uniform(cos_theta_min, cos_theta_max);
+    Double_t theta    = std::acos(costheta);
+    Double_t px       = p * std::cos(phi) * std::sin(theta);
+    Double_t py       = p * std::sin(phi) * std::sin(theta);
+    Double_t pz       = p * std::cos(theta);
+
     // Generates random vectors, uniformly distributed over the surface of a
     // sphere of given radius, in this case momentum.
-    r1->Sphere(px, py, pz, p);
+    //r1->Sphere(px, py, pz, p);
 
     // type 1 is final state
     // pdgid 111 - pi0 135 MeV/c^2
     GenParticlePtr p3 = std::make_shared<GenParticle>(
         FourVector(
             px, py, pz,
-            sqrt((px * px) + (py * py) + (pz * pz) + (0.134976*0.134976))),
+            sqrt(p*p + (0.134976*0.134976))),
         111, 1);
 
     GenVertexPtr v1 = std::make_shared<GenVertex>();
diff --git a/emcal_pi0_reader.cxx b/emcal_pi0_reader.cxx
index d22ce6e4b20444c7136bb6596967470ff7119286..53ddc1776670d7e7c310ce30ffdf1d41efea4d7e 100644
--- a/emcal_pi0_reader.cxx
+++ b/emcal_pi0_reader.cxx
@@ -25,9 +25,9 @@ void emcal_pi0_reader(){
   	gStyle->SetPadGridX(1);
   	gStyle->SetPadGridY(1);
   	gStyle->SetPadLeftMargin(0.14);
-	gStyle->SetPadRightMargin(0.14);
+	gStyle->SetPadRightMargin(0.17);
 
-	ReaderAscii hepmc_input("./data/emcal_pi0_0GeVto30GeV_100kEvt.hepmc");
+	ReaderAscii hepmc_input("./data/emcal_pi0_0GeVto30GeV_1000kEvt.hepmc");
 	int         events_parsed = 0;
 	GenEvent    evt(Units::GEV, Units::MM);
 
@@ -39,7 +39,8 @@ void emcal_pi0_reader(){
 	TH2F* h_pi0_pzpt   = new TH2F("h_pi0_pzpt",    "pi0 pt vs pz;pt [GeV];pz [GeV]",    100,-0.5,30.5,100,-30.5,30.5);
 	TH2F* h_pi0_pxpy   = new TH2F("h_pi0_pxpy",    "pi0 px vs py;px [GeV];py [GeV]",    100,-30.5,30.5,100,-30.5,30.5);
 	TH3F* h_pi0_p      = new TH3F("h_pi0_p",       "pi0 p;px [GeV];py [GeV];pz [GeV]",  100,-30.5,30.5,100,-30.5,30.5,100,-30.5,30.5);
-
+	TH1F* h_pi0_mom    = new TH1F("h_pi0_mom",     "pi0 momentum;p [GeV];Events",       100,-0.5,30.5);
+	
 	while(!hepmc_input.failed()) {
 		// Read event from input file
 		hepmc_input.read_event(evt);
@@ -56,6 +57,9 @@ void emcal_pi0_reader(){
 				h_pi0_pzpt->Fill(TMath::Sqrt(p->momentum().px()*p->momentum().px() + p->momentum().py()*p->momentum().py()),p->momentum().pz());
 				h_pi0_pxpy->Fill(p->momentum().px(),p->momentum().py());
 				h_pi0_p->Fill(p->momentum().px(),p->momentum().py(),p->momentum().pz());
+				h_pi0_mom->Fill(TMath::Sqrt(p->momentum().px()*p->momentum().px() + 
+							p->momentum().py()*p->momentum().py() + 
+							p->momentum().pz()*p->momentum().pz()));
 				}
       			}
     		}
@@ -91,7 +95,7 @@ void emcal_pi0_reader(){
         TCanvas* c3 = new TCanvas("c3","c3", 500, 500);
 	h_pi0_phi->GetYaxis()->SetTitleOffset(1.8);
 	h_pi0_phi->SetLineWidth(2);
-	h_pi0_phi->GetYaxis()->SetRangeUser(0.0,h_pi0_phi->GetMaximum()+100.0);
+	h_pi0_phi->GetYaxis()->SetRangeUser(0.0,h_pi0_phi->GetMaximum()+300.0);
 	h_pi0_phi->SetLineColor(kBlue);
         h_pi0_phi->DrawClone();
 	c3->SaveAs("results/emcal_pi0_phi_0GeVto30GeV.png");
@@ -123,4 +127,12 @@ void emcal_pi0_reader(){
 	c6->SaveAs("results/emcal_pi0_p_0GeVto30GeV.png");
 	c6->SaveAs("results/emcal_pi0_p_0GeVto30GeV.pdf");
 
+        TCanvas* c7 = new TCanvas("c7","c7", 500, 500);
+        h_pi0_mom->GetYaxis()->SetTitleOffset(1.8);
+        h_pi0_mom->SetLineWidth(2);
+        h_pi0_mom->SetLineColor(kBlue);
+        h_pi0_mom->DrawClone();
+        c7->SaveAs("results/emcal_pi0_mom_0GeVto30GeV.png");
+        c7->SaveAs("results/emcal_pi0_mom_0GeVto30GeV.pdf");
+
 }