From 112eb22594dcc8d29e03f7499aa367ee882674f4 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 21 Dec 2018 23:56:38 +0000
Subject: [PATCH] added rotation hack
---
Processes/Sibyll/Interaction.h | 39 ++++++++++++++++++++++++++--------
1 file changed, 30 insertions(+), 9 deletions(-)
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index fb1334da2..23af2bf2a 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -128,11 +128,22 @@ namespace corsika::process::sibyll {
<< process::sibyll::CanInteract(p.GetPID()) << endl;
if (process::sibyll::CanInteract(p.GetPID())) {
const CoordinateSystem& rootCS =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
Point pOrig = p.GetPosition();
TimeType tOrig = p.GetTime();
-
+ // 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);
+
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
@@ -143,7 +154,7 @@ namespace corsika::process::sibyll {
*/
// FOR NOW: set target to oxygen
const int kTarget = corsika::particles::Oxygen::
- GetNucleusA(); // env.GetTargetParticle().GetPID();
+ GetNucleusA(); // env.GetTargetParticle().GetPID();
// FOR NOW: target is always at rest
const hep::MassType nucleon_mass = 0.5*(corsika::particles::Proton::GetMass()+corsika::particles::Neutron::GetMass());
@@ -152,6 +163,9 @@ namespace corsika::process::sibyll {
cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV
<< endl;
+ cout << "beam momentum in sibyll frame (GeV/c): " << p.GetMomentum().GetComponents(sibyllCS) / 1_GeV
+ << endl;
+
cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "time: " << tOrig << endl;
@@ -210,15 +224,15 @@ namespace corsika::process::sibyll {
// add particles from sibyll to stack
// link to sibyll stack
+ // here we need to pass projectile momentum and energy to define the local sibyll frame
+ // and the boosts to the lab. frame
SibStack ss;
- // SibStack does not know about momentum yet so we need counter to access
// momentum array in Sibyll
- int i = -1;
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
EnergyType E_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
- ++i;
+
// skip particles that have decayed in Sibyll
if( psib.HasDecayed()) continue;
@@ -226,7 +240,15 @@ namespace corsika::process::sibyll {
// compute beta_vec * p_vec
// arbitrary Lorentz transformation based on sibyll routines
const auto gammaBetaComponents = gambet.GetComponents();
- const auto pSibyllComponents = psib.GetMomentum().GetComponents();
+ // FOR NOW: fill vector in sibCS and then rotate into rootCS
+ // can be done in SibStack by passing sibCS
+ // get momentum vector in sibyllCS
+ const auto pSibyllComponentsSibCS = psib.GetMomentum().GetComponents();
+ // temporary vector in sibyllCS
+ auto SibVector = MomentumVector( sibyllCS, pSibyllComponentsSibCS);
+ // rotatate to rootCS
+ const auto pSibyllComponents = SibVector.GetComponents(rootCS);
+ // boost to lab. frame
hep::EnergyType en_lab = 0. * 1_GeV;
hep::MomentumType p_lab_components[3];
en_lab = psib.GetEnergy() * gamma;
@@ -245,7 +267,6 @@ namespace corsika::process::sibyll {
auto pnew = s.NewParticle();
pnew.SetEnergy(en_lab);
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
-
corsika::geometry::QuantityVector<energy_hep_d> p_lab_c{
p_lab_components[0], p_lab_components[1], p_lab_components[2]};
pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
--
GitLab