diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index db52594fc2a13c3df511711847264aa318b2a532..2877990e99d0d8a9a120ca3f792e8ea50c66f85d 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -226,12 +226,22 @@ int main() {
   // setup particle stack, and add primary particle
   setup::Stack stack;
   stack.Clear();
-  const hep::EnergyType E0 = 100_GeV;
+  const hep::EnergyType E0 = 1000_GeV;
+  double theta = 45.;
+  double phi = 20.;
   {
     auto particle = stack.NewParticle();
     particle.SetPID(Code::Proton);
     hep::MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV);
-    auto plab = stack::super_stupid::MomentumVector(rootCS, 0_GeV, 0_GeV, -P0);
+    auto momentumComponents = [](double theta, double phi, MomentumType&ptot)
+			      {
+				return std::make_tuple( ptot*sin(theta)*cos(phi), ptot*sin(theta)*sin(phi), ptot*cos(theta) );
+			      };
+    auto const [px, py, pz] = momentumComponents( theta / 180.* M_PI, phi / 180.* M_PI, P0);
+      //    auto plab = stack::super_stupid::MomentumVector(rootCS, 0_GeV, 0_GeV, -P0);
+    auto plab = stack::super_stupid::MomentumVector(rootCS, {px, py, pz});
+    cout << "input angles: theta=" << theta << " phi=" << phi << endl;
+    cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
     particle.SetEnergy(E0);
     particle.SetMomentum(plab);
     particle.SetTime(0_ns);
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 23af2bf2a3ebf4f919c5f8b063a37e2017c1370c..9f1062115b1df4d66a6a53864f1577c5c566f6fc 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -135,14 +135,23 @@ namespace corsika::process::sibyll {
 	// sibyll CS has z along particle momentum
 	// FOR NOW: hard coded z-axis for corsika frame
 	QuantityVector<length_d> const zAxis{0_m, 0_m, 1_m};
-	QuantityVector<length_d> const xAxis{1_m, 0_m, 0_m};
-	auto pt = []( MomentumVector p ){
-		    return sqrt(p.GetComponents()[0] * p.GetComponents()[0] + p.GetComponents()[1] * p.GetComponents()[1]);
-		  };
-	double theta = acos( p.GetMomentum().GetComponents()[2] / p.GetMomentum().norm());
-	cout << "ProcessSibyll: polar angle between sibyllCS and rootCS: " << theta << endl;
-	double phi = asin( p.GetMomentum().GetComponents()[0]/pt(p.GetMomentum() ) );
-	CoordinateSystem sibyllCS = rootCS.rotate(xAxis, theta);
+	QuantityVector<length_d> const yAxis{0_m, 1_m, 0_m};
+	auto rotation_angles = [](MomentumVector const &pin)
+			       {
+				 const auto p = pin.GetComponents();
+				 const auto th =  acos( p[2] / p.norm() );
+				 const auto ph =  atan2( p[1]/1_GeV, p[0]/1_GeV );//acos( p[0] / sqrt(p[0]*p[0]+p[1]*p[1] ) ); 
+				   return std::make_tuple(th, ph);
+			       };
+	// auto pt = []( MomentumVector &p ){
+	// 	    return sqrt(p.GetComponents()[0] * p.GetComponents()[0] + p.GetComponents()[1] * p.GetComponents()[1]);
+	// 	  };
+	// double theta = acos( p.GetMomentum().GetComponents()[2] / p.GetMomentum().norm());
+	auto const [theta, phi] = rotation_angles( p.GetMomentum() );
+	cout << "ProcessSibyll: zenith angle between sibyllCS and rootCS: " << theta / M_PI * 180. << endl;
+	cout << "ProcessSibyll: azimuth angle between sibyllCS and rootCS: " << phi / M_PI * 180. << endl;
+	//double phi = asin( p.GetMomentum().GetComponents()[0]/pt(p.GetMomentum() ) );
+	CoordinateSystem sibyllCS = rootCS.rotate(zAxis,phi).rotate(yAxis, theta);
 	
         /*
            the target should be defined by the Environment,