From 5845823f00c55887e46642ab98a6c8a8efd2028b Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 26 Dec 2018 16:46:55 +0000
Subject: [PATCH] added injection of primary particle with zentith and azimuth
angle, fixed rotation in azimuth
---
Documentation/Examples/cascade_example.cc | 14 +++++++++++--
Processes/Sibyll/Interaction.h | 25 +++++++++++++++--------
2 files changed, 29 insertions(+), 10 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index db52594f..2877990e 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 23af2bf2..9f106211 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,
--
GitLab