diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f1e437ef1a4c84ba42ead78f7b5b509aa2b401c8..fef09923e475538032c86fbac341119799dcd3a9 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -208,13 +208,15 @@ int main() {
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
+ // fraction of oxygen
+ const float fox = 0.20946;
using MyHomogeneousModel =
corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
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,corsika::particles::Code::Proton},
- std::vector<float>{0.5,0.3,0.2}));
+ 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));
@@ -237,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.;
{
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index ff5477f75d5eb5ba9906944059160c352e100098..8ed9ca78d2a8d64f590a6dbcf2918da459487605 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -30,7 +30,6 @@ namespace corsika::process::sibyll {
int fCount = 0;
int fNucCount = 0;
- int kTarget = 0;
public:
Interaction(corsika::environment::Environment const& env) : fEnvironment(env) { }
~Interaction() {
@@ -71,54 +70,6 @@ namespace corsika::process::sibyll {
// read from cross section code table
const int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
-
- // get target from environment
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
- const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
- const auto sample_target = [](const corsika::environment::NuclearComposition& comp)
- {
- double prev =0.;
- std::vector<float> cumFrac;
- for (auto f: comp.GetFractions()){
- cumFrac.push_back(prev+f);
- prev = cumFrac.back();
- }
- int a = 1;
- const double r = s_rndm_(a);
- int i = -1;
- //cout << "rndm: " << r << endl;
- for (auto cf: cumFrac){
- ++i;
- //cout << "cf: " << cf << " " << comp.GetComponents()[i] << endl;
- if(r<cf) break;
- }
- return comp.GetComponents()[i];
- };
-
- const auto compo = currentNode->GetModelProperties().GetNuclearComposition();
- for (auto tid: compo.GetComponents()){
- cout << "Interaction: possible targets: " << tid << endl;
- }
- const auto tg = sample_target( compo );
- cout << "Interaction: target selected: " << tg << endl;
- // test if valid in sibyll
- // FOR NOW: throw error if not, may need workaround for atmosphere, contains small argon fraction
- if(IsNucleus(tg))
- kTarget = GetNucleusA(tg);
- else
- if(tg==corsika::particles::Proton::GetCode())
- kTarget = 1;
- else
- throw std::runtime_error("Sibyll target particle unknown. Only nuclei with A<18 or protons are allowed.");
-
- cout << "Interaction: kTarget: " << kTarget << endl;
- // redundant check
- if(kTarget>18||kTarget==0)
- throw std::runtime_error("Sibyll target outside range. Only nuclei with A<18 are allowed.");
const hep::MassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
corsika::particles::Neutron::GetMass());
@@ -139,33 +90,70 @@ namespace corsika::process::sibyll {
<< " input energy: " << p.GetEnergy() / 1_GeV << endl
<< " beam can interact:" << kBeam << endl
<< " beam XS code:" << kBeam << endl
- << " beam pid:" << p.GetPID() << endl
- << " target mass number:" << kTarget << std::endl;
+ << " beam pid:" << p.GetPID() << endl;
if (kInteraction) {
- double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
- double dumdif[3];
-
- if (kTarget == 1)
- sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
- else
- sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy);
-
- std::cout << "Interaction: "
- << "MinStep: sibyll return: " << prodCrossSection << std::endl;
- si::CrossSectionType sig = prodCrossSection * 1_mbarn;
- std::cout << "Interaction: "
- << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
-
- const si::MassType nucleon_mass =
- 0.93827_GeV / corsika::units::constants::cSquared;
- std::cout << "Interaction: "
- << "nucleon mass " << nucleon_mass << std::endl;
- // calculate interaction length in medium
- GrammageType int_length = kTarget * nucleon_mass / sig;
- std::cout << "Interaction: "
- << "interaction length (g/cm2): " << int_length << std::endl;
+ // get target from environment
+ /*
+ the target should be defined by the Environment,
+ ideally as full particle object so that the four momenta
+ and the boosts can be defined..
+ */
+ const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ const auto mediumComposition = currentNode->GetModelProperties().GetNuclearComposition();
+ // determine average interaction length
+ // FOR NOW: weighted sum
+ int i =-1;
+ double avgTargetMassNumber = 0.;
+ si::CrossSectionType weightedProdCrossSection = 0_mbarn;
+ int kTarget = 0;
+ // get weights of components from environment/medium
+ const auto w = mediumComposition.GetFractions();
+ // loop over components in medium
+ for (auto targetId: mediumComposition.GetComponents() ){
+ i++;
+ cout << "Interaction: get interaction lenght for target: " << targetId << endl;
+ if(IsNucleus(targetId))
+ kTarget = GetNucleusA(targetId);
+ else
+ if(targetId==corsika::particles::Proton::GetCode())
+ kTarget = 1;
+ else
+ throw std::runtime_error("GetInteractionLength: Sibyll target particle unknown. Only nuclei or protons are allowed.");
+
+ cout << "Interaction: kTarget: " << kTarget << endl;
+ // check if nuclei in range
+ if(kTarget>18||kTarget==0)
+ throw std::runtime_error("Sibyll target outside range. Only nuclei with A<18 are allowed.");
+
+ // get cross section from sibyll
+ double sigProd, dummy, dum1, dum2, dum3, dum4;
+ double dumdif[3];
+
+ if (kTarget == 1)
+ sib_sigma_hp_(kBeam, Ecm, dum1, dum2, sigProd, dumdif, dum3, dum4);
+ else
+ sib_sigma_hnuc_(kBeam, kTarget, Ecm, sigProd, dummy);
+
+ std::cout << "Interaction: "
+ << " IntLength: sibyll return: " << sigProd << std::endl;
+ weightedProdCrossSection += w[i] * sigProd * 1_mbarn;
+ avgTargetMassNumber += w[i] * kTarget;
+ }
+ cout << "Interaction: "
+ << "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mbarn
+ << endl;
+
+ const si::MassType nucleon_mass =
+ 0.93827_GeV / corsika::units::constants::cSquared;
+ std::cout << "Interaction: "
+ << "nucleon mass " << nucleon_mass << std::endl;
+ // calculate interaction length in medium
+#warning check interaction length units
+ GrammageType int_length = avgTargetMassNumber * nucleon_mass / weightedProdCrossSection;
+ std::cout << "Interaction: "
+ << "interaction length (g/cm2): " << int_length / ( 0.001_kg ) * 1_cm * 1_cm << std::endl;
return int_length;
}
@@ -219,10 +207,41 @@ namespace corsika::process::sibyll {
const CoordinateSystem tempCS = rootCS.rotate(zAxis, phi);
const CoordinateSystem sibyllCS = tempCS.rotate(yAxis, theta);
-
- // redundant check if target code is valid
- if(kTarget>18||kTarget==0)
- throw std::runtime_error("Interaction: Sibyll target outside range. Only nuclei with A<18 are allowed.");
+ const auto sample_target = [](const corsika::environment::NuclearComposition& comp)
+ {
+ double prev =0.;
+ std::vector<float> cumFrac;
+ for (auto f: comp.GetFractions()){
+ cumFrac.push_back(prev+f);
+ prev = cumFrac.back();
+ }
+ int a = 1;
+ const double r = s_rndm_(a);
+ int i = -1;
+ //cout << "rndm: " << r << endl;
+ for (auto cf: cumFrac){
+ ++i;
+ //cout << "cf: " << cf << " " << comp.GetComponents()[i] << endl;
+ if(r<cf) break;
+ }
+ return comp.GetComponents()[i];
+ };
+ const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ const auto mediumComposition = currentNode->GetModelProperties().GetNuclearComposition();
+
+ const auto tg = sample_target( mediumComposition );
+ cout << "Interaction: target selected: " << tg << endl;
+ // test if valid in sibyll
+ // FOR NOW: throw error if not, may need workaround for atmosphere, contains small argon fraction
+ int kTarget = 0;
+ if(IsNucleus(tg))
+ kTarget = GetNucleusA(tg);
+ else
+ if(tg==corsika::particles::Proton::GetCode())
+ kTarget = 1;
+ else
+ throw std::runtime_error("Sibyll target particle unknown. Only nuclei with A<18 or protons are allowed.");
+
// FOR NOW: target is always at rest
const hep::MassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +