From f8e82c66d4c0a7031024981ac3868db6dc5702ec Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 21 Dec 2018 16:01:50 +0000
Subject: [PATCH] use nucleon mass for center-of-mass energy calculation
---
Documentation/Examples/cascade_example.cc | 2 +-
Processes/Sibyll/Interaction.h | 7 ++++---
2 files changed, 5 insertions(+), 4 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 232e7d00..db52594f 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -226,7 +226,7 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const hep::EnergyType E0 = 1_TeV;
+ const hep::EnergyType E0 = 100_GeV;
{
auto particle = stack.NewParticle();
particle.SetPID(Code::Proton);
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index eb45a570..8fd27958 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -64,8 +64,8 @@ namespace corsika::process::sibyll {
// FOR NOW: assume target is oxygen
const int kTarget = corsika::particles::Oxygen::GetNucleusA();
- hep::EnergyType Etot =
- p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
+ const hep::MassType nucleon_mass = 0.5*(corsika::particles::Proton::GetMass()+corsika::particles::Neutron::GetMass());
+ hep::EnergyType Etot = p.GetEnergy() + nucleon_mass;
MomentumVector Ptot(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
@@ -146,7 +146,8 @@ namespace corsika::process::sibyll {
GetNucleusA(); // env.GetTargetParticle().GetPID();
// FOR NOW: target is always at rest
- const EnergyType Etarget = 0_GeV + corsika::particles::Proton::GetMass();
+ const hep::MassType nucleon_mass = 0.5*(corsika::particles::Proton::GetMass()+corsika::particles::Neutron::GetMass());
+ const EnergyType Etarget = 0_GeV + nucleon_mass;
const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV
--
GitLab