moved target sampling from sibyll process to medium

Documentation/Examples/cascade_example.cc

@@ -239,7 +239,7 @@ int main() {
   // setup particle stack, and add primary particle
   setup::Stack stack;
   stack.Clear();
-  const hep::EnergyType E0 = 100_GeV;
+  const hep::EnergyType E0 = 100_TeV;
   double theta = 0.;
   double phi = 0.;
   {
@@ -269,6 +269,8 @@ int main() {
 
   cout << "Result: E0=" << E0 / 1_GeV << endl;
   cut.ShowResults();
+  const hep::EnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
   cout << "total energy (GeV): "
-       << (cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy()) / 1_GeV << endl;
+       << Efinal / 1_GeV << endl
+       << "relative difference (%): " << (Efinal / E0 - 1. ) * 100 << endl;
 }

Environment/CMakeLists.txt

@@ -22,6 +22,7 @@ target_link_libraries (
   CORSIKAgeometry
   CORSIKAparticles
   CORSIKAunits
+  CORSIKArandom
   )
 
 target_include_directories (

Environment/HomogeneousMedium.h

@@ -17,6 +17,7 @@
 #include <corsika/geometry/Trajectory.h>
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/random/RNGManager.h>
 
 /**
  * a homogeneous medium
@@ -41,6 +42,29 @@ namespace corsika::environment {
     }
     NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
 
+    corsika::particles::Code const& GetTarget( std::vector<corsika::units::si::CrossSectionType> &sigma) const override {
+      using namespace corsika::units::si;
+      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 );
+      }
+					 
+      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() );
+      static corsika::random::RNG& kRNG =
+	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 {
@@ -52,7 +76,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

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

Processes/Sibyll/Interaction.h

@@ -248,37 +248,19 @@ namespace corsika::process::sibyll {
 	// sample target mass number
 	const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
 	const auto mediumComposition = currentNode->GetModelProperties().GetNuclearComposition();
-	// get cross sections and nucleon number for target materials	
-	std::vector<si::CrossSectionType> cross_section_components;
+	// get cross sections for target materials
+	/*
+	  Here we read the cross section from the interaction model again, 
+	  should be passed from GetInteractionLength if possible
+	 */
+#warning reading interaction cross section again, should not be necessary
+	std::vector<si::CrossSectionType> cross_section_of_components;
 	for(auto targetId: mediumComposition.GetComponents() ){
 	  const auto [sigProd, nNuc] = GetCrossSection( corsikaBeamId, targetId, Ecm);	  
-	  cross_section_components.push_back( sigProd );
+	  cross_section_of_components.push_back( sigProd );
 	}	
 
-	// this routine could be moved to the environment as GetTarget( std::vector<si::CrossSectionType> )
-	const auto sample_target = [](const corsika::environment::NuclearComposition& comp, const std::vector<si::CrossSectionType> & sigma_comp )
-				       {
-					 int i=-1;
-					 si::CrossSectionType total_weighted_sigma = 0._mbarn;
-					 std::vector<float> fractions;				 
-					 for(auto w: comp.GetFractions() ){
-					   i++;
-					   cout << "fraction: " << w << endl;
-					   total_weighted_sigma +=  w * sigma_comp[i];
-					   fractions.push_back( w * sigma_comp[i] / 1_mbarn );
-					 }
-					 
-					 for(auto f: fractions){
-					   f = f / ( total_weighted_sigma / 1_mbarn );
-					   cout << "reweighted fraction: " << f << endl;
-					 }
-					 std::discrete_distribution channelDist( fractions.begin(), fractions.end() );
-					 static corsika::random::RNG& kRNG =
-					   corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
-					 const int ichannel = channelDist(kRNG);
-					 return comp.GetComponents()[ichannel];
-				       };
-	const auto targetCode = sample_target( mediumComposition, cross_section_components );
+	const auto targetCode = currentNode->GetModelProperties().GetTarget( cross_section_of_components);
 	cout << "Interaction: target selected: " << targetCode << endl;
 	/*
 	  FOR NOW: allow nuclei with A<18 or protons only.