IAP GITLAB

Skip to content
Snippets Groups Projects
Commit ed38f9de authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

Merge branch 'barn' into 'master'

correct barn literals: _Xbarn -> _Xb

See merge request !97
parents 4fd1a81f f3ca4d59
No related branches found
No related tags found
1 merge request!97correct barn literals
Pipeline #415 passed
...@@ -158,7 +158,7 @@ namespace phys { ...@@ -158,7 +158,7 @@ namespace phys {
QUANTITY_DEFINE_SCALING_LITERALS(eV, hepenergy_d, 1) QUANTITY_DEFINE_SCALING_LITERALS(eV, hepenergy_d, 1)
QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::sigma_d, QUANTITY_DEFINE_SCALING_LITERALS(b, corsika::units::si::sigma_d,
magnitude(corsika::units::constants::barn)) magnitude(corsika::units::constants::barn))
} // namespace literals } // namespace literals
......
...@@ -59,7 +59,7 @@ TEST_CASE("PhysicalUnits", "[Units]") { ...@@ -59,7 +59,7 @@ TEST_CASE("PhysicalUnits", "[Units]") {
REQUIRE(1_mol / 1_amol == Approx(1e18)); REQUIRE(1_mol / 1_amol == Approx(1e18));
REQUIRE(1_K / 1_zK == Approx(1e21)); REQUIRE(1_K / 1_zK == Approx(1e21));
REQUIRE(1_K / 1_yK == Approx(1e24)); REQUIRE(1_K / 1_yK == Approx(1e24));
REQUIRE(1_barn / 1_mbarn == Approx(1e3)); REQUIRE(1_b / 1_mb == Approx(1e3));
REQUIRE(1_A / 1_hA == Approx(1e-2)); REQUIRE(1_A / 1_hA == Approx(1e-2));
REQUIRE(1_m / 1_km == Approx(1e-3)); REQUIRE(1_m / 1_km == Approx(1e-3));
......
...@@ -56,7 +56,7 @@ namespace corsika::process::HadronicElasticModel { ...@@ -56,7 +56,7 @@ namespace corsika::process::HadronicElasticModel {
auto const projectileEnergy = p.GetEnergy(); auto const projectileEnergy = p.GetEnergy();
auto const avgCrossSection = [&]() { auto const avgCrossSection = [&]() {
CrossSectionType avgCrossSection = 0_barn; CrossSectionType avgCrossSection = 0_b;
for (size_t i = 0; i < fractions.size(); ++i) { for (size_t i = 0; i < fractions.size(); ++i) {
auto const targetMass = particles::GetMass(components[i]); auto const targetMass = particles::GetMass(components[i]);
...@@ -65,7 +65,7 @@ namespace corsika::process::HadronicElasticModel { ...@@ -65,7 +65,7 @@ namespace corsika::process::HadronicElasticModel {
avgCrossSection += CrossSection(s) * fractions[i]; avgCrossSection += CrossSection(s) * fractions[i];
} }
std::cout << "avgCrossSection: " << avgCrossSection / 1_mbarn << " mb" std::cout << "avgCrossSection: " << avgCrossSection / 1_mb << " mb"
<< std::endl; << std::endl;
return avgCrossSection; return avgCrossSection;
...@@ -200,8 +200,8 @@ namespace corsika::process::HadronicElasticModel { ...@@ -200,8 +200,8 @@ namespace corsika::process::HadronicElasticModel {
units::si::detail::static_pow<2>(sigmaTotal) / units::si::detail::static_pow<2>(sigmaTotal) /
(16 * M_PI * ConvertHEPToSI<CrossSectionType::dimension_type>(B(s))); (16 * M_PI * ConvertHEPToSI<CrossSectionType::dimension_type>(B(s)));
std::cout << "HEM sigmaTot = " << sigmaTotal / 1_mbarn << " mb" << std::endl; std::cout << "HEM sigmaTot = " << sigmaTotal / 1_mb << " mb" << std::endl;
std::cout << "HEM sigmaElastic = " << sigmaElastic / 1_mbarn << " mb" << std::endl; std::cout << "HEM sigmaElastic = " << sigmaElastic / 1_mb << " mb" << std::endl;
return sigmaElastic; return sigmaElastic;
} }
......
...@@ -156,14 +156,14 @@ namespace corsika::process::pythia { ...@@ -156,14 +156,14 @@ namespace corsika::process::pythia {
const double sigEla = fSigma.sigmaEl(); const double sigEla = fSigma.sigmaEl();
const double sigProd = fSigma.sigmaTot() - sigEla; const double sigProd = fSigma.sigmaTot() - sigEla;
return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn); return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
} else } else
throw std::runtime_error("pythia cross section init failed"); throw std::runtime_error("pythia cross section init failed");
} else { } else {
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn, return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
std::numeric_limits<double>::infinity() * 1_mbarn); std::numeric_limits<double>::infinity() * 1_mb);
} }
} else { } else {
throw std::runtime_error("invalid target for pythia"); throw std::runtime_error("invalid target for pythia");
...@@ -221,7 +221,7 @@ namespace corsika::process::pythia { ...@@ -221,7 +221,7 @@ namespace corsika::process::pythia {
// determine average interaction length // determine average interaction length
// weighted sum // weighted sum
int i = -1; int i = -1;
si::CrossSectionType weightedProdCrossSection = 0_mbarn; si::CrossSectionType weightedProdCrossSection = 0_mb;
// get weights of components from environment/medium // get weights of components from environment/medium
const auto w = mediumComposition.GetFractions(); const auto w = mediumComposition.GetFractions();
// loop over components in medium // loop over components in medium
...@@ -235,13 +235,13 @@ namespace corsika::process::pythia { ...@@ -235,13 +235,13 @@ namespace corsika::process::pythia {
elaCrossSection; // ONLY TO AVOID COMPILER WARNING elaCrossSection; // ONLY TO AVOID COMPILER WARNING
cout << "Interaction: IntLength: pythia return (mb): " cout << "Interaction: IntLength: pythia return (mb): "
<< productionCrossSection / 1_mbarn << endl << productionCrossSection / 1_mb << endl
<< "Interaction: IntLength: weight : " << w[i] << endl; << "Interaction: IntLength: weight : " << w[i] << endl;
weightedProdCrossSection += w[i] * productionCrossSection; weightedProdCrossSection += w[i] * productionCrossSection;
} }
cout << "Interaction: IntLength: weighted CrossSection (mb): " cout << "Interaction: IntLength: weighted CrossSection (mb): "
<< weightedProdCrossSection / 1_mbarn << endl << weightedProdCrossSection / 1_mb << endl
<< "Interaction: IntLength: average mass number: " << "Interaction: IntLength: average mass number: "
<< mediumComposition.GetAverageMassNumber() << endl; << mediumComposition.GetAverageMassNumber() << endl;
......
...@@ -114,7 +114,7 @@ namespace corsika::process::sibyll { ...@@ -114,7 +114,7 @@ namespace corsika::process::sibyll {
sigProd = std::numeric_limits<double>::infinity(); sigProd = std::numeric_limits<double>::infinity();
sigEla = std::numeric_limits<double>::infinity(); sigEla = std::numeric_limits<double>::infinity();
} }
return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn); return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
} }
template <> template <>
...@@ -169,7 +169,7 @@ namespace corsika::process::sibyll { ...@@ -169,7 +169,7 @@ namespace corsika::process::sibyll {
// determine average interaction length // determine average interaction length
// weighted sum // weighted sum
int i = -1; int i = -1;
si::CrossSectionType weightedProdCrossSection = 0_mbarn; si::CrossSectionType weightedProdCrossSection = 0_mb;
// get weights of components from environment/medium // get weights of components from environment/medium
const auto& w = mediumComposition.GetFractions(); const auto& w = mediumComposition.GetFractions();
// loop over components in medium // loop over components in medium
...@@ -183,13 +183,13 @@ namespace corsika::process::sibyll { ...@@ -183,13 +183,13 @@ namespace corsika::process::sibyll {
elaCrossSection; // ONLY TO AVOID COMPILER WARNING elaCrossSection; // ONLY TO AVOID COMPILER WARNING
cout << "Interaction: " cout << "Interaction: "
<< " IntLength: sibyll return (mb): " << productionCrossSection / 1_mbarn << " IntLength: sibyll return (mb): " << productionCrossSection / 1_mb
<< endl; << endl;
weightedProdCrossSection += w[i] * productionCrossSection; weightedProdCrossSection += w[i] * productionCrossSection;
} }
cout << "Interaction: " cout << "Interaction: "
<< "IntLength: weighted CrossSection (mb): " << "IntLength: weighted CrossSection (mb): "
<< weightedProdCrossSection / 1_mbarn << endl; << weightedProdCrossSection / 1_mb << endl;
// calculate interaction length in medium // calculate interaction length in medium
//#warning check interaction length units //#warning check interaction length units
......
...@@ -113,8 +113,8 @@ namespace corsika::process::sibyll { ...@@ -113,8 +113,8 @@ namespace corsika::process::sibyll {
auto const protonId = Code::Proton; auto const protonId = Code::Proton;
auto const [siginel, sigela] = auto const [siginel, sigela] =
fHadronicInteraction.GetCrossSection(protonId, protonId, Ecm); fHadronicInteraction.GetCrossSection(protonId, protonId, Ecm);
const double dsig = siginel / 1_mbarn; const double dsig = siginel / 1_mb;
const double dsigela = sigela / 1_mbarn; const double dsigela = sigela / 1_mb;
// loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile // loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
for (int j = 1; j < fMaxNucleusAProjectile; ++j) { for (int j = 1; j < fMaxNucleusAProjectile; ++j) {
const int jj = j + 1; const int jj = j + 1;
...@@ -169,7 +169,7 @@ namespace corsika::process::sibyll { ...@@ -169,7 +169,7 @@ namespace corsika::process::sibyll {
cout << "ReadCrossSectionTable: " << ia << " " << ib << " " << e0 << endl; cout << "ReadCrossSectionTable: " << ia << " " << ib << " " << e0 << endl;
signuc2_(ia, ib, e0, sig); signuc2_(ia, ib, e0, sig);
cout << "ReadCrossSectionTable: sig=" << sig << endl; cout << "ReadCrossSectionTable: sig=" << sig << endl;
return sig * 1_mbarn; return sig * 1_mb;
} }
// TODO: remove elastic cross section? // TODO: remove elastic cross section?
...@@ -201,13 +201,13 @@ namespace corsika::process::sibyll { ...@@ -201,13 +201,13 @@ namespace corsika::process::sibyll {
if (fHadronicInteraction.IsValidTarget(TargetId)) { if (fHadronicInteraction.IsValidTarget(TargetId)) {
auto const sigProd = ReadCrossSectionTable(iBeamA, TargetId, LabEnergyPerNuc); auto const sigProd = ReadCrossSectionTable(iBeamA, TargetId, LabEnergyPerNuc);
cout << "cross section (mb): " << sigProd / 1_mbarn << endl; cout << "cross section (mb): " << sigProd / 1_mb << endl;
return std::make_tuple(sigProd, 0_mbarn); return std::make_tuple(sigProd, 0_mb);
} else { } else {
throw std::runtime_error("target outside range."); throw std::runtime_error("target outside range.");
} }
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn, return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
std::numeric_limits<double>::infinity() * 1_mbarn); std::numeric_limits<double>::infinity() * 1_mb);
} }
template <> template <>
...@@ -277,7 +277,7 @@ namespace corsika::process::sibyll { ...@@ -277,7 +277,7 @@ namespace corsika::process::sibyll {
// determine average interaction length // determine average interaction length
// weighted sum // weighted sum
int i = -1; int i = -1;
si::CrossSectionType weightedProdCrossSection = 0_mbarn; si::CrossSectionType weightedProdCrossSection = 0_mb;
// get weights of components from environment/medium // get weights of components from environment/medium
const auto w = mediumComposition.GetFractions(); const auto w = mediumComposition.GetFractions();
// loop over components in medium // loop over components in medium
...@@ -290,13 +290,13 @@ namespace corsika::process::sibyll { ...@@ -290,13 +290,13 @@ namespace corsika::process::sibyll {
[[maybe_unused]] auto elaCrossSectionCopy = elaCrossSection; [[maybe_unused]] auto elaCrossSectionCopy = elaCrossSection;
cout << "NuclearInteraction: " cout << "NuclearInteraction: "
<< "IntLength: nuclib return (mb): " << productionCrossSection / 1_mbarn << "IntLength: nuclib return (mb): " << productionCrossSection / 1_mb
<< endl; << endl;
weightedProdCrossSection += w[i] * productionCrossSection; weightedProdCrossSection += w[i] * productionCrossSection;
} }
cout << "NuclearInteraction: " cout << "NuclearInteraction: "
<< "IntLength: weighted CrossSection (mb): " << "IntLength: weighted CrossSection (mb): "
<< weightedProdCrossSection / 1_mbarn << endl; << weightedProdCrossSection / 1_mb << endl;
// calculate interaction length in medium // calculate interaction length in medium
GrammageType const int_length = mediumComposition.GetAverageMassNumber() * GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
...@@ -472,8 +472,8 @@ namespace corsika::process::sibyll { ...@@ -472,8 +472,8 @@ namespace corsika::process::sibyll {
const auto protonId = particles::Proton::GetCode(); const auto protonId = particles::Proton::GetCode();
const auto [prodCrossSection, elaCrossSection] = const auto [prodCrossSection, elaCrossSection] =
fHadronicInteraction.GetCrossSection(protonId, protonId, EcmNN); fHadronicInteraction.GetCrossSection(protonId, protonId, EcmNN);
const double sigProd = prodCrossSection / 1_mbarn; const double sigProd = prodCrossSection / 1_mb;
const double sigEla = elaCrossSection / 1_mbarn; const double sigEla = elaCrossSection / 1_mb;
// sample number of interactions (only input variables, output in common cnucms) // sample number of interactions (only input variables, output in common cnucms)
// nuclear multiple scattering according to glauber (r.i.p.) // nuclear multiple scattering according to glauber (r.i.p.)
int_nuc_(kATarget, kAProj, sigProd, sigEla); int_nuc_(kATarget, kAProj, sigProd, sigEla);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment