From 0278fd67042bec7dadfb9ec04b76ece295841c61 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Fri, 18 Jan 2019 19:19:16 +0100
Subject: [PATCH] cascade
---
Documentation/Examples/CMakeLists.txt | 3 ++-
Documentation/Examples/cascade_example.cc | 9 ++++++-
.../HadronicElasticModel.h | 27 ++++++++++++-------
3 files changed, 28 insertions(+), 11 deletions(-)
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 00de097e7..30073cfdc 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -25,10 +25,11 @@ add_executable (cascade_example cascade_example.cc)
target_compile_options(cascade_example PRIVATE -g) # do not skip asserts
target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogging
CORSIKArandom
- ProcessSibyll
+ #~ ProcessSibyll
CORSIKAcascade
ProcessStackInspector
ProcessTrackWriter
+ ProcessHadronicElasticModel
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 9e44b146b..bb1dfd37b 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -12,6 +12,7 @@
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -235,10 +236,15 @@ int main() {
corsika::process::sibyll::Interaction sibyll(env);
corsika::process::sibyll::Decay decay;
ProcessCut cut(8_GeV);
+
+ corsika::random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
+ corsika::process::HadronicElasticModel::HadronicElasticInteraction hadronicElastic(env);
+
corsika::process::TrackWriter::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list
- auto sequence = p0 << sibyll << decay << cut << trackWriter;
+ //~ auto sequence = p0 << sibyll << decay << cut << trackWriter;
+ auto sequence = hadronicElastic << trackWriter;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
@@ -246,6 +252,7 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
+
const HEPEnergyType E0 = 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
double theta = 0.;
double phi = 0.;
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h
index e96cc7703..243d9727e 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.h
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.h
@@ -21,6 +21,8 @@
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
+#include <iomanip>
+
namespace corsika::process::HadronicElasticModel {
/**
* A simple model for elastic hadronic interactions based on the formulas
@@ -164,7 +166,7 @@ namespace corsika::process::HadronicElasticModel {
auto const sqrtS = eProjectileCoM + eTargetCoM;
auto const s = units::si::detail::static_pow<2>(sqrtS);
-
+
auto const B = this->B(sqrtS);
std::cout << B << std::endl;
@@ -189,19 +191,26 @@ namespace corsika::process::HadronicElasticModel {
<< ')' << std::endl;
auto const theta =
- 2 *
- asin(sqrt(absT / (4 * pProjectileCoMSqNorm))); // would work for lab frame, too
+ 2 * asin(sqrt(absT / (4 * pProjectileCoMSqNorm))); // would work for in any frame
auto const phi = phiDist(fRNG);
geometry::QuantityVector<units::si::hepmomentum_d> const scatteredMomentum{
pProjectileCoMNorm * sin(theta) * cos(phi),
- pProjectileCoMNorm * sin(theta) * cos(phi), pProjectileCoMNorm * cos(theta)};
- auto const [eProjectileScattered, pProjectileScattered] =
+ pProjectileCoMNorm * sin(theta) * sin(phi), pProjectileCoMNorm * cos(theta)};
+
+ auto const [eProjectileScatteredLab, pProjectileScatteredLab] =
boost.fromCoM(eProjectileCoM, scatteredMomentum);
-
- p.SetMomentum(pProjectileScattered);
- p.SetEnergy(sqrt(pProjectileScattered.squaredNorm() +
- units::si::detail::static_pow<2>(targetMass)));
+
+ std::cout << "xxx: " << pProjectileScatteredLab.squaredNorm() / p.GetMomentum().squaredNorm() << std::endl;
+
+ p.SetMomentum(pProjectileScatteredLab);
+
+ // alternatives: 1. do not touch energy
+ // 2. take back-boosted energy
+ // 3. recalculate energy from momentum
+
+ //~ p.SetEnergy(sqrt(pProjectileScatteredLab.squaredNorm() +
+ //~ units::si::detail::static_pow<2>(corsika::particles::GetMass(p.GetPID()))));
return process::EProcessReturn::eOk;
}
--
GitLab