Environment/NuclearComposition.h

@@ -16,6 +16,7 @@
 #include <corsika/units/PhysicalUnits.h>
 
 #include <cassert>
+#include <functional>
 #include <numeric>
 #include <random>
 #include <stdexcept>
@@ -82,7 +83,6 @@ namespace corsika::environment {
     auto WeightedSum(TFunction func) const {
       using ResultQuantity = decltype(func(*fComponents.cbegin()));
 
-      auto const sum = [](auto x, auto y) { return x + y; };
       auto const prod = [&](auto const compID, auto const fraction) {
         return func(compID) * fraction;
       };
@@ -91,12 +91,12 @@ namespace corsika::environment {
         return std::inner_product(
             fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
             ResultQuantity::zero(), // .zero() is defined for quantity types only
-            sum, prod);
+            std::plus<ResultQuantity>(), prod);
       } else {
         return std::inner_product(
             fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
             ResultQuantity(0), // in other cases we have to use a bare 0
-            sum, prod);
+            std::plus<ResultQuantity>(), prod);
       }
     }
 

Processes/Pythia/Interaction.cc

@@ -76,37 +76,37 @@ namespace corsika::process::pythia {
   }
 
   void Interaction::SetParticleListStable(
-      std::vector<particles::Code> const& particleList) {
-    for (auto p : particleList) Interaction::SetStable(p);
+      std::vector<particles::Code> const& vParticleList) {
+    for (auto p : vParticleList) Interaction::SetStable(p);
   }
 
-  void Interaction::SetUnstable(const particles::Code pCode) {
-    cout << "Pythia::Interaction: setting " << pCode << " unstable.." << endl;
-    fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), true);
+  void Interaction::SetUnstable(const particles::Code vCode) {
+    cout << "Pythia::Interaction: setting " << vCode << " unstable.." << endl;
+    fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(vCode)), true);
   }
 
-  void Interaction::SetStable(const particles::Code pCode) {
-    cout << "Pythia::Interaction: setting " << pCode << " stable.." << endl;
-    fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), false);
+  void Interaction::SetStable(const particles::Code vCode) {
+    cout << "Pythia::Interaction: setting " << vCode << " stable.." << endl;
+    fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(vCode)), false);
   }
 
   void Interaction::ConfigureLabFrameCollision(
-      const particles::Code BeamId, const particles::Code TargetId,
-      const units::si::HEPEnergyType BeamEnergy) {
+      const particles::Code vBeamId, const particles::Code vTargetId,
+      const units::si::HEPEnergyType vBeamEnergy) {
     using namespace units::si;
     // Pythia configuration of the current event
     // very clumsy. I am sure this can be done better..
 
     // set beam
     // beam id for pythia
-    auto const pdgBeam = static_cast<int>(particles::GetPDG(BeamId));
+    auto const pdgBeam = static_cast<int>(particles::GetPDG(vBeamId));
     std::stringstream stBeam;
     stBeam << "Beams:idA = " << pdgBeam;
     fPythia.readString(stBeam.str());
     // set target
-    auto pdgTarget = static_cast<int>(particles::GetPDG(TargetId));
+    auto pdgTarget = static_cast<int>(particles::GetPDG(vTargetId));
     // replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
-    if (TargetId == particles::Code::Hydrogen)
+    if (vTargetId == particles::Code::Hydrogen)
       pdgTarget = static_cast<int>(particles::GetPDG(particles::Code::Proton));
     std::stringstream stTarget;
     stTarget << "Beams:idB = " << pdgTarget;
@@ -114,7 +114,7 @@ namespace corsika::process::pythia {
     // set frame to lab. frame
     fPythia.readString("Beams:frameType = 2");
     // set beam energy
-    const double Elab = BeamEnergy / 1_GeV;
+    const double Elab = vBeamEnergy / 1_GeV;
     std::stringstream stEnergy;
     stEnergy << "Beams:eA = " << Elab;
     fPythia.readString(stEnergy.str());
@@ -124,28 +124,28 @@ namespace corsika::process::pythia {
     fPythia.init();
   }
 
-  bool Interaction::CanInteract(const corsika::particles::Code pCode) {
-    return pCode == corsika::particles::Code::Proton ||
-           pCode == corsika::particles::Code::Neutron ||
-           pCode == corsika::particles::Code::AntiProton ||
-           pCode == corsika::particles::Code::AntiNeutron ||
-           pCode == corsika::particles::Code::PiMinus ||
-           pCode == corsika::particles::Code::PiPlus;
+  bool Interaction::CanInteract(const corsika::particles::Code vCode) {
+    return vCode == corsika::particles::Code::Proton ||
+           vCode == corsika::particles::Code::Neutron ||
+           vCode == corsika::particles::Code::AntiProton ||
+           vCode == corsika::particles::Code::AntiNeutron ||
+           vCode == corsika::particles::Code::PiMinus ||
+           vCode == corsika::particles::Code::PiPlus;
   }
 
   tuple<units::si::CrossSectionType, units::si::CrossSectionType>
-  Interaction::GetCrossSection(const particles::Code BeamId,
-                               const particles::Code TargetId,
-                               const units::si::HEPEnergyType CoMenergy) {
+  Interaction::GetCrossSection(const particles::Code vBeamId,
+                               const particles::Code vTargetId,
+                               const units::si::HEPEnergyType vCoMenergy) {
     using namespace units::si;
 
     // interaction possible in pythia?
     if (TargetId == particles::Code::Proton || TargetId == particles::Code::Hydrogen) {
-      if (CanInteract(BeamId) && ValidCoMEnergy(CoMenergy)) {
+      if (CanInteract(vBeamId) && ValidCoMEnergy(vCoMenergy)) {
         // input particle PDG
-        auto const pdgCodeBeam = static_cast<int>(particles::GetPDG(BeamId));
-        auto const pdgCodeTarget = static_cast<int>(particles::GetPDG(TargetId));
-        const double ecm = CoMenergy / 1_GeV;
+        auto const pdgCodeBeam = static_cast<int>(particles::GetPDG(vBeamId));
+        auto const pdgCodeTarget = static_cast<int>(particles::GetPDG(vTargetId));
+        const double ecm = vCoMenergy / 1_GeV;
 
         // calculate cross section
         fSigma.calc(pdgCodeBeam, pdgCodeTarget, ecm);
@@ -168,7 +168,7 @@ namespace corsika::process::pythia {
   }
 
   template <>
-  units::si::GrammageType Interaction::GetInteractionLength(Particle& p, Track&) {
+  units::si::GrammageType Interaction::GetInteractionLength(Particle& vP, Track&) {
 
     using namespace units;
     using namespace units::si;
@@ -178,7 +178,7 @@ namespace corsika::process::pythia {
     CoordinateSystem& rootCS =
         RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
-    const particles::Code corsikaBeamId = p.GetPID();
+    const particles::Code corsikaBeamId = vP.GetPID();
 
     // beam particles for pythia : 1, 2, 3 for p, pi, k
     // read from cross section code table
@@ -188,9 +188,9 @@ namespace corsika::process::pythia {
     process::pythia::MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
 
     // total momentum and energy
-    HEPEnergyType Elab = p.GetEnergy() + constants::nucleonMass;
+    HEPEnergyType Elab = vP.GetEnergy() + constants::nucleonMass;
     process::pythia::MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
-    pTotLab += p.GetMomentum();
+    pTotLab += vP.GetMomentum();
     pTotLab += pTarget;
     auto const pTotLabNorm = pTotLab.norm();
     // calculate cm. energy
@@ -198,9 +198,9 @@ namespace corsika::process::pythia {
         (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
 
     cout << "Interaction: LambdaInt: \n"
-         << " input energy: " << p.GetEnergy() / 1_GeV << endl
+         << " input energy: " << vP.GetEnergy() / 1_GeV << endl
          << " beam can interact:" << kInteraction << endl
-         << " beam pid:" << p.GetPID() << endl;
+         << " beam pid:" << vP.GetPID() << endl;
 
     // TODO: move limits into variables
     if (kInteraction && Elab >= 8.5_GeV && ValidCoMEnergy(ECoM)) {
@@ -211,30 +211,16 @@ namespace corsika::process::pythia {
         ideally as full particle object so that the four momenta
         and the boosts can be defined..
       */
-      const auto* currentNode = p.GetNode();
+      const auto* currentNode = vP.GetNode();
       const auto mediumComposition =
           currentNode->GetModelProperties().GetNuclearComposition();
       // determine average interaction length
-      // weighted sum
-      int i = -1;
-      si::CrossSectionType weightedProdCrossSection = 0_mb;
-      // get weights of components from environment/medium
-      const auto& w = mediumComposition.GetFractions();
-      // loop over components in medium
-      for (auto const targetId : mediumComposition.GetComponents()) {
-        i++;
-        cout << "Interaction: get interaction length for target: " << targetId << endl;
-
-        auto const [productionCrossSection, elaCrossSection] =
-            GetCrossSection(corsikaBeamId, targetId, ECoM);
-        [[maybe_unused]] const auto& dummy_elaCrossSection = elaCrossSection;
-
-        cout << "Interaction: IntLength: pythia return (mb): "
-             << productionCrossSection / 1_mb << endl
-             << "Interaction: IntLength: weight : " << w[i] << endl;
-
-        weightedProdCrossSection += w[i] * productionCrossSection;
-      }
+
+      auto const weightedProdCrossSection =
+          mediumComposition.WeightedSum([=](auto vTargetID) {
+            return std::get<0>(this->GetCrossSection(corsikaBeamId, vTargetID, ECoM));
+          });
+
       cout << "Interaction: IntLength: weighted CrossSection (mb): "
            << weightedProdCrossSection / 1_mb << endl
            << "Interaction: IntLength: average mass number: "

Processes/Pythia/Interaction.h

@@ -37,9 +37,9 @@ namespace corsika::process::pythia {
     void SetStable(const corsika::particles::Code);
 
     bool WasInitialized() { return fInitialized; }
-    bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
+    bool ValidCoMEnergy(corsika::units::si::HEPEnergyType vEcm) {
       using namespace corsika::units::si;
-      return (10_GeV < ecm) && (ecm < 1_PeV);
+      return (10_GeV < vEcm) && (vEcm < 1_PeV);
     }
 
     bool CanInteract(const corsika::particles::Code);
@@ -47,9 +47,9 @@ namespace corsika::process::pythia {
                                     const corsika::particles::Code,
                                     const corsika::units::si::HEPEnergyType);
     std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
-    GetCrossSection(const corsika::particles::Code BeamId,
-                    const corsika::particles::Code TargetId,
-                    const corsika::units::si::HEPEnergyType CoMenergy);
+    GetCrossSection(const corsika::particles::Code,
+                    const corsika::particles::Code,
+                    const corsika::units::si::HEPEnergyType);
 
     template <typename TParticle, typename TTrack>
     corsika::units::si::GrammageType GetInteractionLength(TParticle&, TTrack&);