xMerge branch 'sibyll' of gitlab.ikp.kit.edu:AirShowerPhysics/corsika into sibyll

Documentation/Examples/cascade_example.cc

@@ -70,7 +70,7 @@ public:
     sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
     
     std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
-    CrossSectionType sig = prodCrossSection  / 1000. * barn;
+    CrossSectionType sig = prodCrossSection * 1_mbarn;
     std::cout << "ProcessSplit: " << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
 
     const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;

Framework/Units/PhysicalConstants.h

@@ -54,8 +54,9 @@ namespace corsika::units::si::constants {
   // unified atomic mass unit
   constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
 
-  // barn
-  constexpr quantity<area_d> barn{Rep(1.e-28L) * meter * meter};
+  // barn moved to PhysicalUnits
+  //  constexpr quantity<area_d> barn{Rep(1.e-28L) * meter * meter};
+  
   
   // etc.
 

Framework/Units/PhysicalUnits.h

@@ -24,17 +24,6 @@ namespace phys {
   }   // namespace units
 } // namespace phys
 
-namespace phys {
-  namespace units {
-    namespace literals {
-      QUANTITY_DEFINE_SCALING_LITERALS(barn, area_d,
-                                       magnitude(corsika::units::si::constants::barn))
-
-      // phys::units::quantity<energy_d/mass_d> Joule2Kg = c2; // 1_Joule / 1_kg;
-
-    } // namespace literals
-  }   // namespace units
-}
 
 namespace phys {
   namespace units {
@@ -49,6 +38,7 @@ namespace phys {
 }
 
 
+
 namespace corsika::units::si {
   using namespace phys::units;
   using namespace phys::units::literals;
@@ -63,10 +53,45 @@ namespace corsika::units::si {
   using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
   using MassType = phys::units::quantity<phys::units::mass_d, double>;
 
-  using CrossSectionType = phys::units::quantity<phys::units::area_d, double>;
+  // defining momentum you suckers
+  // dimensions, i.e. composition in base SI dimensions
+  using momentum_d                     = phys::units::dimensions< 1, 1, -1 >;
+  // defining the unit of momentum, so far newton-meter, maybe go to HEP?
+  constexpr phys::units::quantity< momentum_d > newton_second   { meter * kilogram / second };
+  // defining the type
+  using MomentumType = phys::units::quantity<momentum_d, double>;
+
+  // defining cross section
+  using sigma_d      = phys::units::dimensions< 2, 0, 0 >;
+  constexpr phys::units::quantity< sigma_d  > barn       {Rep(1.e-28L) * meter * meter};
+  using CrossSectionType = phys::units::quantity<sigma_d, double>;
   
 } // end namespace corsika::units::si
 
+namespace phys {
+  namespace units {
+    namespace literals {
+      QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::sigma_d,
+                                       magnitude(corsika::units::si::constants::barn))
+
+      // phys::units::quantity<energy_d/mass_d> Joule2Kg = c2; // 1_Joule / 1_kg;
+
+    } // namespace literals
+  }   // namespace units
+}
+
+namespace phys {
+  namespace units {
+    namespace literals {
+      QUANTITY_DEFINE_SCALING_LITERALS(newton_second, corsika::units::si::momentum_d,
+                                       magnitude(corsika::units::si::newton_second))
+
+      // phys::units::quantity<energy_d/mass_d> Joule2Kg = c2; // 1_Joule / 1_kg;
+
+    } // namespace literals
+  }   // namespace units
+}
+
 // we want to call the operator<< without namespace... I think
 using namespace phys::units::io;
 

Framework/Units/testUnits.cc

@@ -32,7 +32,7 @@ TEST_CASE("PhysicalUnits", "[Units]") {
 
     LengthType l1 = 10_nm;
     l1 = l1;
-
+    
     LengthType arr0[5];
     arr0[0] = 5_m;
 
@@ -41,6 +41,9 @@ TEST_CASE("PhysicalUnits", "[Units]") {
     std::array<EnergyType, 4> arr2; // empty array
 
     [[maybe_unused]] std::array<EnergyType, 4> arr3 = {1_GeV, 1_eV, 5_MeV};
+
+    [[maybe_unused]] auto p1 = 10_newton_second;
+    REQUIRE(p1 == 10_newton_second);
   }
 
   SECTION("Powers in literal units") {
@@ -54,7 +57,7 @@ TEST_CASE("PhysicalUnits", "[Units]") {
     REQUIRE(1_mol / 1_amol == Approx(1e18));
     REQUIRE(1_K / 1_zK == Approx(1e21));
     REQUIRE(1_K / 1_yK == Approx(1e24));
-    //    REQUIRE(1_barn / 1_mbarn == Approx(1e3));
+    REQUIRE(1_barn / 1_mbarn == Approx(1e3));
     
     REQUIRE(1_A / 1_hA == Approx(1e-2));
     REQUIRE(1_m / 1_km == Approx(1e-3));
@@ -77,6 +80,10 @@ TEST_CASE("PhysicalUnits", "[Units]") {
     REQUIRE(E2 == 40_GeV);
     REQUIRE(E2 / 1_GeV == Approx(40));
 
+    const MassType m = 1_kg;
+    const SpeedType v = 1_m / 1_s;
+    REQUIRE( m*v == 1_newton_second);
+    
     const double lgE = log10(E2 / 1_GeV);
     REQUIRE(lgE == Approx(log10(40.)));