IAP GITLAB

Skip to content
Snippets Groups Projects
Commit c6e4bdba authored by Felix Riehn's avatar Felix Riehn
Browse files

extended cross section app for multi target-proj combinations, hack job you can do better than this

parent f787b4b2
No related branches found
No related tags found
1 merge request!649Resolve "hadronic interactions with Pythia"
Pipeline #14212 passed
Stage: quality
Stage: config
Stage: build_test
Stage: build_test_example
Stage: optional
Hide whitespace changes
Inline Side-by-side
examples/physics_examples/cross_section.cpp
+ 50
− 58
View file @ c6e4bdba
...@@ -44,22 +44,24 @@ using namespace corsika; ...@@ -44,22 +44,24 @@ using namespace corsika;
* *
* @param model the interaction model to use * @param model the interaction model to use
* @param model_name string with model name * @param model_name string with model name
* @param lab_frame switch between lab. frame or center-of-mass input
* *
* @return a tuple of: inelastic cross section, elastic cross section * @return a tuple of: inelastic cross section, elastic cross section
*/ */
template <typename TModel> template <typename TModel>
void calculate_cross_sections(TModel& model, std::string const& model_name, void calculate_cross_sections(TModel& model, std::string const& model_name) {
bool lab_frame = true) {
std::stringstream fname; std::stringstream fname;
fname << "cross-sections-" << model_name << ".txt"; fname << "cross-sections-" << model_name << ".txt";
std::ofstream out(fname.str()); std::ofstream out(fname.str());
out << "# cross section table for " << model_name << "\n"; out << "# cross section table for " << model_name << "\n";
out << "# i, E_lab (GeV), E_cm per nucleon (GeV), cross section p-p (mb), production " out << "# i, P_lab (GeV), E_cm per nucleon (GeV), cross section p-p (mb), cross "
"section pi-p (mb), production "
"cross section p-O16 (mb), production cross section p-N14 (mb), production " "cross section p-O16 (mb), production cross section p-N14 (mb), production "
"cross section p-Ar40 (mb) \n"; "cross section p-Ar40 (mb), cross section pi-O16 (mb), production cross "
"section, pi-N14 (mb), production, cross section pi-Ar40 (mb), cross section "
"He-O16 (mb), production cross section, He-N14 (mb), production, cross section "
"He-Ar40 (mb)\n";
CoordinateSystemPtr const& cs = get_root_CoordinateSystem(); CoordinateSystemPtr const& cs = get_root_CoordinateSystem();
...@@ -71,62 +73,51 @@ void calculate_cross_sections(TModel& model, std::string const& model_name, ...@@ -71,62 +73,51 @@ void calculate_cross_sections(TModel& model, std::string const& model_name,
for (int i = 0; i < nebins; ++i) { for (int i = 0; i < nebins; ++i) {
corsika::units::si::HEPMomentumType const setP = pow(10, da * i + amin) * 1_GeV; corsika::units::si::HEPMomentumType const setP = pow(10, da * i + amin) * 1_GeV;
CrossSectionType xs_prod_hNuc[3][3] = {};
CrossSectionType xs_prod_hp[3] = {};
HEPEnergyType setE, comEnn;
int l = -1;
for (auto projId : {Code::Proton, Code::PiPlus, Code::Helium}) {
++l;
setE = calculate_total_energy(setP, get_mass(projId));
comEnn = corsika::calculate_com_energy(setE, get_mass(projId),
corsika::constants::nucleonMass);
auto setE = calculate_total_energy(setP, Proton::mass);
corsika::units::si::HEPEnergyType const comEnn = corsika::calculate_com_energy(
setE, Proton::mass, corsika::constants::nucleonMass);
FourMomentum p4protonProj(0_eV, {cs, 0_GeV, 0_GeV, 0_GeV});
FourMomentum p4protonTarg(0_eV, {cs, 0_GeV, 0_GeV, 0_GeV});
FourMomentum p4oxygenTarg(0_eV, {cs, 0_GeV, 0_GeV, 0_GeV});
FourMomentum p4nitrogenTarg(0_eV, {cs, 0_GeV, 0_GeV, 0_GeV});
FourMomentum p4argonTarg(0_eV, {cs, 0_GeV, 0_GeV, 0_GeV});
if (lab_frame) {
// four momenta in lab frame // four momenta in lab frame
p4protonProj = corsika::FourMomentum(setE, MomentumVector(cs, {0_eV, 0_eV, setP})); corsika::FourMomentum p4Proj(setE, MomentumVector(cs, {0_eV, 0_eV, setP}));
p4protonTarg = corsika::FourMomentum p4protonTarg(Proton::mass,
corsika::FourMomentum(Proton::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV})); MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
p4oxygenTarg =
corsika::FourMomentum(Oxygen::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV})); // had-p
p4nitrogenTarg = auto const [xs_prod, xs_ela] =
corsika::FourMomentum(Nitrogen::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV})); model->getCrossSectionInelEla(projId, Code::Proton, p4Proj, p4protonTarg);
p4argonTarg = xs_prod_hp[l] = xs_prod;
corsika::FourMomentum(Argon::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
} else { int k = -1;
// four momenta in com frame for (auto targNuc : {Code::Oxygen, Code::Nitrogen, Code::Argon}) {
p4protonProj = ++k;
corsika::FourMomentum(corsika::calculate_total_energy(setP / 2, Proton::mass), FourMomentum const p4nucTarg = corsika::FourMomentum(
MomentumVector(cs, {0_eV, 0_eV, setP / 2})); get_mass(targNuc), MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
p4protonTarg = xs_prod_hNuc[l][k] = model->getCrossSection(projId, targNuc, p4Proj, p4nucTarg);
corsika::FourMomentum(corsika::calculate_total_energy(setP / 2, Proton::mass), }
MomentumVector(cs, {0_eV, 0_eV, -setP / 2}));
p4oxygenTarg =
corsika::FourMomentum(corsika::calculate_total_energy(setP / 2, Oxygen::mass),
MomentumVector(cs, {0_eV, 0_eV, -setP / 2}));
} }
// p-p
auto const [xs_prod_pp, xs_ela_pp] = model->getCrossSectionInelEla(
Code::Proton, Code::Proton, p4protonProj, p4protonTarg);
// p-O
CrossSectionType const xs_prod_pO =
model->getCrossSection(Code::Proton, Code::Oxygen, p4protonProj, p4oxygenTarg);
CrossSectionType const xs_prod_pN = model->getCrossSection(
Code::Proton, Code::Nitrogen, p4protonProj, p4nitrogenTarg);
CrossSectionType const xs_prod_pAr =
model->getCrossSection(Code::Proton, Code::Argon, p4protonProj, p4argonTarg);
CORSIKA_LOG_INFO( CORSIKA_LOG_INFO(
"Elab={:15.2f} GeV,Ecom={:15.2f} GeV, sig-p-p={:6.2f} mb, sig-p-O={:6.2f} mb, " "Elab={:15.2f} GeV,Ecom={:15.2f} GeV, sig-p-p={:6.2f} mb, sig-pi-p={:6.2f} mb, "
"sig-p-N={:6.2f} mb, sig-p-Ar={:6.2f} mb", "p-O={:6.2f} mb, p-N={:6.2f} mb, p-Ar={:6.2f} mb "
setE / 1_GeV, comEnn / 1_GeV, xs_prod_pp / 1_mb, xs_prod_pO / 1_mb, "pi-O={:6.2f} mb, pi-N={:6.2f} mb, pi-Ar={:6.2f} mb "
xs_prod_pN / 1_mb, xs_prod_pAr / 1_mb); "He-O={:6.2f} mb, He-N={:6.2f} mb, He-Ar={:6.2f} mb",
setP / 1_GeV, comEnn / 1_GeV, xs_prod_hp[0] / 1_mb, xs_prod_hp[1] / 1_mb,
out << i << " " << setE / 1_GeV << " " << comEnn / 1_GeV << " " << xs_prod_pp / 1_mb xs_prod_hNuc[0][0] / 1_mb, xs_prod_hNuc[0][1] / 1_mb, xs_prod_hNuc[0][2] / 1_mb,
<< " " << xs_prod_pO / 1_mb << " " << xs_prod_pN / 1_mb << " " xs_prod_hNuc[1][0] / 1_mb, xs_prod_hNuc[1][1] / 1_mb, xs_prod_hNuc[1][2] / 1_mb,
<< xs_prod_pAr / 1_mb << "\n"; xs_prod_hNuc[2][0] / 1_mb, xs_prod_hNuc[2][1] / 1_mb, xs_prod_hNuc[2][2] / 1_mb);
out << i << " " << setP / 1_GeV << " " << comEnn / 1_GeV << " "
<< xs_prod_hp[0] / 1_mb << " " << xs_prod_hp[1] / 1_mb << " "
<< xs_prod_hNuc[0][0] / 1_mb << " " << xs_prod_hNuc[0][1] / 1_mb << " "
<< xs_prod_hNuc[0][2] / 1_mb << " " << xs_prod_hNuc[1][0] / 1_mb << " "
<< xs_prod_hNuc[1][1] / 1_mb << " " << xs_prod_hNuc[1][2] / 1_mb << " "
<< " " << xs_prod_hNuc[2][0] / 1_mb << " " << xs_prod_hNuc[2][1] / 1_mb << " "
<< xs_prod_hNuc[2][2] / 1_mb << "\n";
} }
out.close(); out.close();
std::cout << "wrote cross section table for " << model_name std::cout << "wrote cross section table for " << model_name
...@@ -148,7 +139,8 @@ int main(int argc, char** argv) { ...@@ -148,7 +139,8 @@ int main(int argc, char** argv) {
if (int_model_name == "sibyll") { if (int_model_name == "sibyll") {
RNGManager<>::getInstance().registerRandomStream("sibyll"); RNGManager<>::getInstance().registerRandomStream("sibyll");
auto model = std::make_shared<corsika::sibyll::InteractionModel>( auto model = std::make_shared<corsika::sibyll::InteractionModel>(
std::set<Code>{Code::Proton, Code::Oxygen}, corsika::setup::C7trackedParticles); std::set<Code>{Code::Proton, Code::Oxygen, Code::Nitrogen},
corsika::setup::C7trackedParticles);
calculate_cross_sections(model, int_model_name); calculate_cross_sections(model, int_model_name);
} else if (int_model_name == "pythia8") { } else if (int_model_name == "pythia8") {
RNGManager<>::getInstance().registerRandomStream("pythia"); RNGManager<>::getInstance().registerRandomStream("pythia");
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum