Merge branch '116-use-central-lorentz-boost-in-sibyll-interface' into 'master'

Resolve "use central Lorentz boost in sibyll interface"

Closes #116

See merge request AirShowerPhysics/corsika!56

Documentation/Examples/cascade_example.cc

@@ -218,8 +218,9 @@ int main() {
   theMedium->SetModelProperties<MyHomogeneousModel>(
       1_kg / (1_m * 1_m * 1_m),
       corsika::environment::NuclearComposition(
-					       std::vector<corsika::particles::Code>{corsika::particles::Code::Nitrogen,corsika::particles::Code::Oxygen},
-					       std::vector<float>{(float)1.-fox, fox}));
+          std::vector<corsika::particles::Code>{corsika::particles::Code::Nitrogen,
+                                                corsika::particles::Code::Oxygen},
+          std::vector<float>{(float)1. - fox, fox}));
 
   universe.AddChild(std::move(theMedium));
 
@@ -272,8 +273,8 @@ int main() {
 
   cout << "Result: E0=" << E0 / 1_GeV << endl;
   cut.ShowResults();
-  const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
-  cout << "total energy (GeV): "
-       << Efinal / 1_GeV << endl
-       << "relative difference (%): " << (Efinal / E0 - 1. ) * 100 << endl;
+  const HEPEnergyType Efinal =
+      cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+  cout << "total energy (GeV): " << Efinal / 1_GeV << endl
+       << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
 }

Environment/HomogeneousMedium.h

@@ -16,8 +16,8 @@
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/Trajectory.h>
 #include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
 #include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
 
 /**
  * a homogeneous medium
@@ -42,29 +42,30 @@ namespace corsika::environment {
     }
     NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
 
-    corsika::particles::Code const& GetTarget( std::vector<corsika::units::si::CrossSectionType> &sigma) const override {
+    corsika::particles::Code const& GetTarget(
+        std::vector<corsika::units::si::CrossSectionType>& sigma) const override {
       using namespace corsika::units::si;
-      int i=-1;
+      int i = -1;
       corsika::units::si::CrossSectionType total_weighted_sigma = 0._mbarn;
-      std::vector<float> fractions;				 
-      for(auto w: fNuclComp.GetFractions() ){
-	i++;
-	std::cout << "HomogeneousMedium: fraction: " << w << std::endl;
-	total_weighted_sigma +=  w * sigma[i];
-	fractions.push_back( w * sigma[i] / 1_mbarn );
+      std::vector<float> fractions;
+      for (auto w : fNuclComp.GetFractions()) {
+        i++;
+        std::cout << "HomogeneousMedium: fraction: " << w << std::endl;
+        total_weighted_sigma += w * sigma[i];
+        fractions.push_back(w * sigma[i] / 1_mbarn);
       }
-					 
-      for(auto f: fractions){
-	f = f / ( total_weighted_sigma / 1_mbarn );
-	std::cout << "HomogeneousMedium: reweighted fraction: " << f << std::endl;
+
+      for (auto f : fractions) {
+        f = f / (total_weighted_sigma / 1_mbarn);
+        std::cout << "HomogeneousMedium: reweighted fraction: " << f << std::endl;
       }
-      std::discrete_distribution channelDist( fractions.begin(), fractions.end() );
+      std::discrete_distribution channelDist(fractions.begin(), fractions.end());
       static corsika::random::RNG& kRNG =
-	corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+          corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
       const int ichannel = channelDist(kRNG);
       return fNuclComp.GetComponents()[ichannel];
     }
-    
+
     corsika::units::si::GrammageType IntegratedGrammage(
         corsika::geometry::Trajectory<corsika::geometry::Line> const&,
         corsika::units::si::LengthType pTo) const override {
@@ -76,7 +77,7 @@ namespace corsika::environment {
         corsika::geometry::Trajectory<corsika::geometry::Line> const&,
         corsika::units::si::GrammageType pGrammage) const override {
       return pGrammage / fDensity;
-    }    
+    }
   };
 
 } // namespace corsika::environment

Environment/IMediumModel.h

@@ -39,7 +39,8 @@ namespace corsika::environment {
 
     virtual NuclearComposition const& GetNuclearComposition() const = 0;
 
-    virtual corsika::particles::Code const& GetTarget( std::vector<corsika::units::si::CrossSectionType> &) const = 0;
+    virtual corsika::particles::Code const& GetTarget(
+        std::vector<corsika::units::si::CrossSectionType>&) const = 0;
   };
 
 } // namespace corsika::environment

Framework/ProcessSequence/ProcessSequence.h

@@ -138,7 +138,8 @@ namespace corsika::process {
 
     template <typename Particle, typename Stack>
     EProcessReturn SelectInteraction(
-        Particle& p, Stack& s, [[maybe_unused]]corsika::units::si::InverseGrammageType lambda_select,
+        Particle& p, Stack& s,
+        [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
         corsika::units::si::InverseGrammageType& lambda_inv_count) {
 
       if constexpr (is_process_sequence<T1type>::value) {
@@ -204,9 +205,10 @@ namespace corsika::process {
 
     // select decay process
     template <typename Particle, typename Stack>
-    EProcessReturn SelectDecay(Particle& p, Stack& s,
-                               [[maybe_unused]] corsika::units::si::InverseTimeType decay_select,
-                               corsika::units::si::InverseTimeType& decay_inv_count) {
+    EProcessReturn SelectDecay(
+        Particle& p, Stack& s,
+        [[maybe_unused]] corsika::units::si::InverseTimeType decay_select,
+        corsika::units::si::InverseTimeType& decay_inv_count) {
       if constexpr (is_process_sequence<T1>::value) {
         // if A is a process sequence --> check inside
         const EProcessReturn ret = A.SelectDecay(p, s, decay_select, decay_inv_count);

Processes/Sibyll/Decay.h

@@ -169,7 +169,7 @@ namespace corsika::process {
         using corsika::geometry::Point;
         using namespace corsika::units::si;
 
-	// TODO: this should be done in a central, common place. Not here..
+        // TODO: this should be done in a central, common place. Not here..
 #ifndef CORSIKA_OSX
         feenableexcept(FE_INVALID);
 #endif
@@ -185,7 +185,8 @@ namespace corsika::process {
         pin.SetMomentum(p.GetMomentum());
         // setting particle mass with Corsika values, may be inconsistent with sibyll
         // internal values
-	// TODO: #warning setting particle mass with Corsika values, may be inconsistent with sibyll internal values
+        // TODO: #warning setting particle mass with Corsika values, may be inconsistent
+        // with sibyll internal values
         pin.SetMass(corsika::particles::GetMass(pCode));
         // remember position
         Point decayPoint = p.GetPosition();
@@ -219,7 +220,7 @@ namespace corsika::process {
         // empty sibyll stack
         ss.Clear();
 
-	// TODO: this should be done in a central, common place. Not here..
+        // TODO: this should be done in a central, common place. Not here..
 #ifndef CORSIKA_OSX
         fedisableexcept(FE_INVALID);
 #endif

Processes/Sibyll/testSibyll.cc

@@ -70,9 +70,9 @@ TEST_CASE("Sibyll", "[processes]") {
 
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/particles/ParticleProperties.h>
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
-#include <corsika/particles/ParticleProperties.h>
 
 #include <corsika/environment/Environment.h>
 #include <corsika/environment/HomogeneousMedium.h>
@@ -83,12 +83,12 @@ using namespace corsika::units;
 
 TEST_CASE("SibyllInterface", "[processes]") {
 
-    // setup environment, geometry
+  // setup environment, geometry
   corsika::environment::Environment env;
   auto& universe = *(env.GetUniverse());
 
   auto theMedium = corsika::environment::Environment::CreateNode<geometry::Sphere>(
-									 geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+      geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
       1_km * std::numeric_limits<double>::infinity());
 
   using MyHomogeneousModel =
@@ -109,14 +109,13 @@ TEST_CASE("SibyllInterface", "[processes]") {
   geometry::Line line(origin, v);
   geometry::Trajectory<geometry::Line> track(line, 10_s);
 
-  
   SECTION("InteractionInterface") {
 
     setup::Stack stack;
     auto particle = stack.NewParticle();
-    
+
     Interaction model(env);
-    
+
     model.Init();
     [[maybe_unused]] const process::EProcessReturn ret =
         model.DoInteraction(particle, stack);
@@ -131,7 +130,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
     {
       const HEPEnergyType E0 = 10_GeV;
       particle.SetPID(particles::Code::Proton);
-      HEPMomentumType P0 = sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
+      HEPMomentumType P0 =
+          sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
       auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
       particle.SetEnergy(E0);
       particle.SetMomentum(plab);
@@ -139,9 +139,9 @@ TEST_CASE("SibyllInterface", "[processes]") {
       geometry::Point p(cs, 0_m, 0_m, 0_m);
       particle.SetPosition(p);
     }
-    
+
     Decay model;
-    
+
     model.Init();
     /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(particle,
                                                                           stack);