cascade

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

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>
@@ -234,10 +235,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";
@@ -245,7 +251,7 @@ int main() {
   // setup particle stack, and add primary particle
   setup::Stack stack;
   stack.Clear();
-  const HEPEnergyType E0 = 100_TeV;
+  const HEPEnergyType E0 = 10_GeV;
   double theta = 0.;
   double phi = 0.;
   {

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;
     }